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Summary 


A 2D linear aeroacoustic theory for rotor/stator interaction with unsteady coupling was 
derived and explored in Volume 1 of this report. Computer program CUP2D has been written 
in FORTRAN embodying the theoretical equations. This volume (Volume 2) describes the 
structure of the code, installation and running, preparation of the input file, and interpretation of 
the output. A sample case is provided with printouts of the input and output. The source code 
is included with comments linking it closely to the theoretical equations in Volume 1. 



Section 1 
Introduction 


This volume provides documentation and user information for the coupled 2D linearized 
cascade code CUP2D. Theory for the code is derived and explored in Volume 1 of this report. 
Material herein discusses how to install and run the code, explains the input file and the printed 
output, outlines the code structure, and provides a listing of the source code. 

CUP2D is written as strictly as possible in FORTRAN 77 and is self-contained so that no 
system subroutines are needed. One exception is that the DOUBLE COMPLEX variable type 
is used; this should be accepted by any modem compiler. The only other exception is that the 
subroutine TIMDAT calls a system-dependent time and date function. This has been found to 

TM TM 

work on Sun and Iris systems but can be deleted or modified by the user, if necessary. 
Figure 1 shows the hierarchy of subroutines with a brief description of the subroutine functions. 
More description can be found in the subroutine comments. To interpret the figure 1, note that 
each routine calls only those routines indented underneath. Thus, for example, READIN calls 
only ALFBET, RTCOEF, GTWAKE, and PRNTIN. Each routine is called only once with the 
exception of GETVS and DSWK, as shown in figure 1. 

The entire code is supplied on disk in a single module (or file) called cup2d.f and can be 
compiled on a UNIX system by entering f77 -o cup2d cup2d.f . This generates an executable 
file which can be run by typing cup2d . The code then looks for the input file cupin.dat , 
which must be in the same directory as cup2d . Normally, the output is written to the screen 
only. However, if the user wishes the output written to a file, cupout.dat for example, the 
command cup2d > cupout.dat would be used. 

Sections 2 and 3 give detailed descriptions of the input and output. The source listing in 
Section 4 is heavily commented with descriptions of subroutine function at the top and 
throughout each routine. Also, to help in linking the code to the theory, variable names were 
chosen to be as close to the names used in the theory derivation as possible. Finally, wherever 
appropriate in the code, equation numbers from Volume 1 are given next to the corresponding 
FORTRAN statements. 
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HIERARCHY OF SUBROUTINES 


PROGRAM CUP2D . . 

READIN 

ALFBET . . . 
RTCOEF . . . 
GTWAKE . . 

PRNTIN 

TIM DAT 


main program 

read input file (UNIT 8) and compute some constants 

compute arrays of alphas and betas 

compute arrays of reflection and transmission coefficients 

compute upwash vectors 

print input 

print time and date of execution 


INFNS 

RCOEFS 

GET VS 

SCOEFS 

GETVS 

GENKRR 

DSWK & WAVE . 

GENKSS 

DSWK & WAVE . 


compute elements of KMATRX 

rotor on stator effect 

compute Smith’s v) , v' 2 , v 3 

stator on rotor effect 

compute Smith’s v) , v' 2 , v 3 

rotor on rotor effect 

Smith’s routines for matrix elements 

stator on stator effect 

Smith’s routines for matrix elements 


SOLVE solve coupled system for loading on both blade rows 

MATINV invert the matrix KMATRX 

LINPAC Routines 

LOADS compute loads from [KMATRX]' 1 * WASH = LOAD 


OUTPUT print out sound pressure and sound power by mode 

GETPWL compute modal sound power 


Figure 1 . This figure indicates all subroutine calls. Except where shown, each routine is 
called only once. 
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Section 2 

Explanation of Input File 


The sample case input for code CUP2D is supplied on the same disk with the source listing 
in the file called cupin.dat , which is shown in figure 2. To facilitate verification of input, all 
of the input numbers are printed with the normal output. Brief definitions of the input are 
included at the bottom of the file. Some further explanations are provided below. 

Line 1 - The comment is provided for user convenience and is printed on the 4 th line of output. 

Line 2 - In the theory, a 2D cascade is considered to be wrapped into a narrow annulus to 
simulate a fan and to permit a mixture of fan nomenclature and cascade nomenclature. In 
particular, this permits numbers of blades to appear directly. The radius to the annulus provides 
a common dimension for the 2 blade rows which is considered to be the effective radius of the 
fan. It is used for non-dimensionalization in the axial spacing of the blade rows. In simulating 
a fan, the effective radius could be taken as 85 percent of the tip radius, in which case the rotor 
rotational Mach numbers at that radius would be input below in line 6. 

Line 3 - The number of panels N p is fixed to be the same for both blade rows. The number 
of harmonics N h is the number used for the coupling equations and in the printout of sound 
pressure and sound power. The code is delivered with dimensioning for a max N h = 5 and a max 
N p = 51. Section 5 gives array dimension information if this needs to be changed. 

Line 4 - See figure 2 of this volume. 

Line 5 - Speed of sound is in feet per second. Density is in pounds mass per cubic foot. 

Line 6 - The code can treat counter-rotation configurations. In this case, input a positive 
rotational Mach number for the front rotor and a negative number for the rear. If either row is 
a stator, set its rotational Mach number to zero. 

Line 7 - Here the user specifies the axial locations where he wants the modal sound pressure to 
be evaluated for the final table in the printout. Distances are measured downstream from the 
front row leading edge and are normalized by rotor effective radius. 

Line 8 and 9 - For line 8, input the value of INTYPE to be used by the the subroutine 
GTWAKE in evaluating the upwash at the two blade rows for excitation of the system. 

INTYPE = 1 is used to apply the Silverstein wake formulas derived in appendix E of 
Volume 1 of this report. Here, the user only specifies the drag coefficient on line 9 and the 
code computes upwash at N control points along the chord for each of N h harmonics. 
This option was used for figure 15 of Volume 1. 

INTYPE = 2 is the same as 1 except that the harmonics above BPF are set to zero. This 
is useful for evaluation of the frequency scattering effect and was used for most of the 
figures in Volume 1. 
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INTYPE = 3 also applies the formulas from appendix E of Volume 1 for computing upwash 
along the stator chord. However, in Equations E-12 and E-14, the absolute values of the 
wake harmonics (F n in those equations) are input directly by the user in line 9. Since the 
wake amplitudes do not decay using this input, this is equivalent to specifying excitation 
by a vorticity wave. INTYPE 3 was used to check the Kousen/Verdon results in Section 
4 of Volume 1. 

INTYPE = 4 is provided so that the user can excite the rotor and stator with an upwash 
distribution of his own choosing. Thus, the upwash vectors WREXT(/i,i) and WSEXT(«,i) 
that would computed from wake formulas using INTYPE 1 are input directly for harmonic 
order n and control point i. These are entered as real numbers in tabular form starting 
on line 9 as shown below for a two harmonic case. 


Real[WREXT( 1,1)] 

Imag[WREXT( 1,1)] 

Real[WREXT( 1 fl p )] 

Imag[WREXT( 1 ,N p )] 

Real[WREXT(2, 1 )] 

Imag[WREXT(2, 1 )] 

Real[WREXT(2 v /V /J )] 

Imag[WREXT(2„/V p )] 

Real[WSEXT(l,l)] 

Imag[WSEXT(l,l)] 

Real [W SEXT( 1 ,N p )] 

Imag[WSEXT(l,A^)] 

Real [ W S EXT (2,1)] 

Imag[WSEXT(2, 1)] 

RealfWSEXT(2,V p )] 

Imag[WSEXT(2,yvp] 


This mode of input could be used to simulate excitation of one blade row by the potential 
field of the other or could be used to simulate blade vibration effects. 
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File cupin.dat 


' Sample case for Code CUP2D, B = 38, V=72' 

38 72 .7628 .9567 .556 

30 3 

.318 .360 .427 0.231 

1037.7 1070.0 1070.3 0.0293 0.0328 0.0326 

.75 0.0 

-0.217 0.86 

2 

.02 

The above is input for a sample case for code CUP2D 

Line 1 Comment - up to 70 characters - in single quotes 

Line 2 # Blades-upst ream row, # Blades-downstream row, 

gap/chord-1, gap/chord-2, 

axial dist LEI to LE2 normalized by fan effective radius 

Line 3 Number of panels each row, Number of harmonics 

Line 4 Axial Mach number - upstream, inter-row, downstream 
Tangential Mach number - inter-row 

Line 5 Speed of sound upstream, inter-row, downstream 
Density upstream, inter-row, downstream 

Line 6 Rotational Mach number- front blade row (0 for IGV) 

Rotational Mach number- rear blade row (0 for EGV, negative # for 
rotor) 

Line 7 Axial locations for acoustic pressure printout normalized by fan 

effective radius, measured positive downstream from front blade leading 
edge . 

Line 8 1 or 2 for input based on drag coefficient. 2 sets the wake harmonics 

above BPF to zero. See code documentation for other options. 

Line 9 Drag coefficient 


Figure 2. Input data set for sample case. 
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Section 3 

Explanation of Code Output 

Output for the sample case is shown in figure 3. Most of the input is printed on the first page. 
Subscripts 1 and 2 refer to the upstream and downstream blade rows respectively. Also, 
subscripts a, b, and c refer to the regions upstream of the upstream row, between rows, and 
downstream of the downstream row. 

For the axial spacing of the blade rows, the user inputs the distance from between leading edges 
of the rows in fan effective radii. The code then computes and prints the axial distance from the 
upstream trailing edge to the downstream leading edge normalized by the upstream chord. Input 
values printed in the output include: 

Mxa, Mxb, and Mxc - axial Mach numbers 
Ms - swirl Mach number (in Region b) 

My], My2 - blade row rotational Mach numbers 
RHOa, RHOb, and RHOc - densities 
Aa, Ab, Ac - speeds of sound 

Relative Mach numbers of the two blade rows Mrell and Mrel2 are computed from the velocity 
triangles in figure 5 of Volume 1. Smith’s reduced frequencies are based on full chord. 

Flow angles Theta 1 and Theta 2 and Swirl Angle are computed from the input Mach numbers. 
Note that Theta 1 is normally negative per figure 5 of Volume 1. 

The long table entitled "EXTERNAL VELOCITY IMPOSED ON CASCADE" gives the 
upwash values used as external excitation of the system. These are listed by harmonic order N 
and control point along the chord, counted by I. The control points are at Smith’s unevenly 
spaced locations given by = 0.5*(J - cos[n(2I-l )/(2NJ]}. Values in the table become the 
vectors WREXT and WSEXT used in the call to the SOLVE routine. 

After the listing of the upwash vectors, the printout in figure 3 shows the items "Entering 
RCOEFS" , and so forth to show the user how near execution is to completion. The "condition 
number" indicates whether the KMATRX is close to being singular. 

"VALUES OF LIFT COEFFICIENTS" are the A C p values computed in the LOADS routine 
integrated over the chords of each blade row. 

In the final table showing modal sound pressures and sound powers, the FREQ column gives the 
value of O nk = nB { M yl - kB 2 M y2 and the mode column gives nBj - kB 2 , which is defined so 
that positive values correspond to co-rotating modes. (Co-rotation implies a mode rotating in the 
direction of positive rotation of the blade rows and in the direction of positive swirl.) The cutoff 
ratios on the right are printed to help with diagnosis. For example, with the first harmonic (N=l) 
the K=1 mode can be seen to be cut on between the rotor and stator (region b) and cut off in the 
upstream and downstream regions. (Cutoff ratios larger than 9.99 are printed as 9.99.) This is 
the "trapped mode" discussed at length in Volume 1. It produces pressure but no power in 
regions a and b. Note that the downstream powers for n = 2 and 3 (namely 56. 1 dB and 62.6 
dB) can be found in the top part of figure 14 in Volume 1 for the rotor rotational Mach number 
= 0.75. 
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Code CUP2D for coupled cascade aeroacoust ics - Version 1.1 

Developed for NASA-Lewis by Pratt & Whitney under Contract NAS3-25952 - Task 10 
Theory documented in NASA CR-4506, Volume I. 

COMMENT: Sample case for Code CUP2D, B^38, V=72 
Time of execution: Mon Mar 1 12:56:37 1993 

Bl= 38 B2= 72 

Gap/Chord(l) = 0.763, Gap/Chord{2) = 0.957 

(Rotor LE to Stator LE) / (Local Radius of Rotor) = 0.556 (input) 

Axial Spacing Between Blade Rows/Rotor Chord 1.9951 (computed) 

Drag Coefficient = 0.020 

Number of panel s= 30, Number of harmonics= 3 

Mxa Mxb Mxc Ms Myl My2 Mr ell Mrel2 

0.318 0.360 0.427 0.231 0.750 0.000 0.632 0.428 

RHOa RHOb RHOc Aa Ab Ac 

0.02930 0.03280 0.03260 1037.7 1070.0 1070.3 

Remainder of printout is computed from input above 
Smiths reduced freqs 0 BPF front row, rear row = 18.532 6.078 

Thetal, Theta2 (in degrees) = -55.253 32.687 

Swirl Angle (in degrees) = 32.69 

Ambient Pressure/Sea Level STD {uptream, downstream) = 0.331 0.391 


EXTERNAL VELOCITY IMPOSED ON CASCADE 


N 

I 

WREXT(real, imag) 

WSEXT ( real , imag) 

1 

1 

0.0000 

0.0000 

0.0186 

0.0261 

1 

2 

0.0000 

0.0000 

0.0194 

0.0254 

1 

3 

0.0000 

0.0000 

0.0210 

0.0241 

1 

4 

0.0000 

0.0000 

0.0232 

0.0218 

1 

5 

0 . 0000 

0.0000 

0.0257 

0.0186 

1 

6 

0.0000 

0 . 0000 

0.0282 

0.0142 

1 

7 

0.0000 

0.0000 

0.0302 

0.0087 

1 

8 

0.0000 

0.0000 

0.0312 

0.0021 

1 

9 

0.0000 

0.0000 

0.0306 

-0.0053 

1 

10 

0.0000 

0.0000 

0.0280 

-0.0128 

1 

11 

0.0000 

0.0000 

0.0233 

-0.0198 

1 

12 

0.0000 

0.0000 

0.0166 

-0 . 0254 

1 

13 

0.0000 

0.0000 

0.0082 

-0 . 0289 

1 

14 

0.0000 

0.0000 

-0.0011 

-0.0298 

1 

15 

0.0000 

0.0000 

-0.0102 

-0.0277 

1 

16 

0.0000 

0.0000 

-0.0182 

-0.0229 

1 

17 

0.0000 

0.0000 

-0.0242 

-0 . 0160 

1 

18 

0.0000 

0.0000 

-0.0277 

-0 . 0078 

1 

19 

0.0000 

0.0000 

-0.0285 

0.0009 

1 

20 

0.0000 

0.0000 

-0.0269 

0 .0089 

1 

21 

0.0000 

0.0000 

-0.0233 

0.0158 

1 

22 

0.0000 

0.0000 

-0 . 0184 

0.0210 

1 

23 

0.0000 

0.0000 

-0.0128 

0.0246 

1 

24 

0.0000 

0.0000 

-0.0073 

0 . 0266 

1 

25 

0.0000 

0.0000 

-0 . 0022 

0.0273 

1 

26 

0.0000 

0.0000 

0.0021 

0.0272 

1 

27 

0.0000 

0.0000 

0.0056 

0.0266 

1 

28 

0.0000 

0.0000 

0.0082 

0.0259 

1 

29 

0.0000 

0.0000 

0.0098 

0.0253 

1 

30 

0.0000 

0.0000 

0.0107 

0.0249 


Figure 3 (beginning). Output for sample case. 
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2 

1 

0.0000 

0.0000 

0.0000 

0.0000 

2 

2 

0.0000 

0.0000 

0.0000 

0 . 0000 

2 

3 

0.0000 

0.0000 

0.0000 

0.0000 

2 

4 

0.0000 

0.0000 

0 . 0000 

0.0000 

2 

5 

0.0000 

0.0000 

0.0000 

0.0000 

2 

6 

0.0000 

0.0000 

0.0000 

0.0000 

2 

7 

0.0000 

0.0000 

0.0000 

0.0000 

2 

8 

0.0000 

0.0000 

0.0000 

0.0000 

2 

9 

0.0000 

0.0000 

0.0000 

0.0000 

2 

10 

0.0000 

0 . 0000 

0.0000 

0 .0000 

2 

11 

0.0000 

0.0000 

0.0000 

0.0000 

2 

12 

0.0000 

0.0000 

0 . 0000 

0.0000 

2 

13 

0.0000 

0.0000 

0.0000 

0.0000 

2 

14 

0.0000 

0 . 0000 

0.0000 

0.0000 

2 

15 

0.0000 

0.0000 

0.0000 

0.0000 

2 

16 

0.0000 

0.0000 

0.0000 

0.0000 

2 

17 

0.0000 

0.0000 

0 . 0000 

0 . 0000 

2 

18 

0.0000 

0.0000 

0 . 0000 

0.0000 

2 

19 

0.0000 

0.0000 

0.0000 

0.0000 

2 

20 

0 . 0000 

0.0000 

0.0000 

0.0000 

2 

21 

0.0000 

0.0000 

0.0000 

0.0000 

2 

22 

0.0000 

0.0000 

0.0000 

0.0000 

2 

23 

0.0000 

0.0000 

0.0000 

0.0000 

2 

24 

0 . 0000 

0.0000 

0.0000 

0.0000 

2 

25 

0.0000 

0 . 0000 

0 . 0000 

0 .0000 

2 

26 

0.0000 

0 . 0000 

0.0000 

0.0000 

2 

27 

0.0000 

0.0000 

0 . 0000 

0 .0000 

2 

28 

0.0000 

0.0000 

0.0000 

0.0000 

2 

29 

0.0000 

0 . 0000 

0.0000 

0 . 0000 

2 

30 

0 . 0000 

0.0000 

0.0000 

0 . 0000 

3 

1 

0 . 0000 

0.0000 

0.0000 

0 .0000 

3 

2 

0.0000 

0.0000 

0.0000 

0 . 0000 

3 

3 

0.0000 

0.0000 

0 . 0000 

0 .0000 

3 

4 

0 . 0000 

0 . 0000 

0.0000 

0 . 0000 

3 

5 

0 . 0000 

0.0000 

0.0000 

0.0000 

3 

6 

0 . 0000 

0.0000 

0.0000 

0.0000 

3 

7 

0.0000 

0.0000 

0 . 0000 

0.0000 

3 

8 

0.0000 

0.0000 

0 . 0000 

0.0000 

3 

9 

0.0000 

0.0000 

0.0000 

0 . 0000 

3 

10 

0.0000 

0.0000 

0 . 0000 

0 . 0000 

3 

11 

0.0000 

0.0000 

0.0000 

0.0000 

3 

12 

0.0000 

0.0000 

0.0000 

0 .0000 

3 

13 

0.0000 

0.0000 

0.0000 

0.0000 

3 

14 

0.0000 

0.0000 

0.0000 

0.0000 

3 

15 

0.0000 

0.0000 

0.0000 

0 . 0000 

3 

16 

0.0000 

0.0000 

0.0000 

0.0000 

3 

17 

0.0000 

0.0000 

0.0000 

0.0000 

3 

18 

0.0000 

0.0000 

0.0000 

0.0000 

3 

19 

0.0000 

0.0000 

0.0000 

0.0000 

3 

20 

0.0000 

0.0000 

0.0000 

0.0000 

3 

21 

0.0000 

0 . 0000 

0.0000 

0.0000 

3 

22 

0.0000 

0 . 0000 

0.0000 

0 . 0000 

3 

23 

0.0000 

0 . 0000 

0 . 0000 

0.0000 

3 

24 

0.0000 

0.0000 

0.0000 

0.0000 

3 

25 

0.0000 

0.0000 

0 .0000 

0.0000 

3 

26 

0.0000 

0.0000 

0.0000 

0.0000 

3 

27 

0.0000 

0.0000 

0.0000 

0.0000 

3 

28 

0.0000 

0.0000 

0.0000 

0.0000 

3 

29 

0.0000 

0.0000 

0.0000 

0.0000 

3 

30 

0.0000 

0.0000 

0.0000 

0.0000 


Figure 3 (continued), output for sample case. 
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Entering RCOEFS 
Entering SCOEFS 
Entering GENKRR 
Entering GENKSS 
Entering MATINV 

Condition number of KMATRX = 4 . 3478710669121D-05 

*** VALUES OF LIFT COEFFICIENTS *** 

N CLROTOR (N) CL STATOR (N) 

real imag real imag 

1 -0.00130 -0.00070 0.00270 -0.01780 

2 0.00009 -0.00004 -0.00029 0.00064 

3 0.00001 0.00003 -0.00032 -0.00051 


Axial locations for sound pressure output in radii from rotor leading edge 
For PRESup, Xa = -.217, For PRESdn, Xc = 0.860 

Decibel Levels for Pressure Waves and Power Levels I Cutoff Ratios 


N 

K 

FREQ 

nBl-kB2 PRESup 

PRESdn 

PWLup 

PWLdn 

1 

A 


B 


c 

1 

-1 

28 . 

SO 

110 

0. 

0. 

0 . 

0 

0. 

0 

0. 

0 

0 . 

,28 

0 . 

,03 

0. 

29 

1 

0 

28. 

50 

38 

0 . 

0 

0 . 

0 

0. 

0 

0. 

0 

0 . 

. 82 

0. 

.56 

0. 

83 

1 

1 

28 . 

,50 

-3 4 

69 . 

9 

21 . 

1 

0. 

.0 

0 . 

.0 

0. 

.91 

1 . 

,15 

0 . 

93 

1 

2 

28. 

.50 

-106 

0. 

0 

0 . 

0 

0. 

,0 

0. 

.0 

0 , 

.29 

0 . 

.54 

0 . 

,30 

— 


— 


Total 

power 

for 

N= 1 


0. 

.0 

0. 

.0 







2 

-3 

57 , 

.00 

292 

0. 

, 0 

0. 

0 

0. 

. 0 

0. 

. 0 

0 

.21 

0 

.04 

0. 

.22 

2 

-2 

57 . 

.00 

220 

0. 

. 0 

0. 

0 

0. 

. 0 

0. 

. 0 

0 

.28 

0 

.03 

0. 

.29 

2 

-1 

57 . 

.00 

148 

0. 

, 0 

0. 

0 

0, 

. 0 

0. 

. 0 

0 

. 42 

0 

.17 

0. 

.43 

2 

0 

57 , 

.00 

76 

0. 

, 0 

0 . 

0 

0. 

. 0 

0. 

. 0 

0 

. 82 

0 

. 56 

0. 

.83 

2 

1 

57 . 

.00 

4 

98. 

.7 

108 . 

9 

41 . 

. 0 

56. 

. 1 

9 

.99 

9 

.99 

9 . 

. 99 

2 

2 

57 , 

.00 

-68 

31 . 

, 6 

0. 

0 

0. 

, 0 

0. 

. 0 

0 

.91 

1 

.15 

0. 

.93 

2 

3 

57 . 

.00 

-140 

0. 

.0 

0 . 

0 

0. 

. 0 

0. 

. 0 

0 

. 44 

0 

. 68 

0. 

. 45 

— 


— 


Total 

power 

for 

N= 2 


41. 

. 0 

56, 

. 1 







3 

-3 

85 

. 50 

330 

0 

.0 

0. 

0 

0. 

.0 

0, 

.0 

0 

.28 

0 

.03 

0 

.29 

3 

-2 

85 

.50 

258 

0 

. 0 

0. 

0 

0, 

. 0 

0, 

.0 

0 

.36 

0 

. 11 

0 

.37 

3 

-1 

85 

.50 

186 

0 

. 0 

0 . 

0 

0, 

.0 

0, 

.0 

0 

. 50 

0 

.25 

0 

.51 

3 

0 

85 

. 50 

114 

0 

. 0 

0. 

0 

0, 

.0 

0 

.0 

0 

.82 

0 

.56 

0 

.83 

3 

1 

85 

. 50 

42 

103 

. 4 

116. 

3 

45 

.4 

62 

.4 

2 

.21 

1 

.93 

2 

.25 

3 

2 

85 

. 50 

-30 

88 

. 9 

102. 

1 

31 

. 1 

48 

.7 

3 

. 10 

3 

.30 

3 

.15 

3 

3 

85 

.50 

-102 

6 

. 5 

0. 

0 

0 

.0 

0 

.0 

0 

. 91 

1 

.15 

0 

.93 

— 


— 


Total 

power 

for 

N= 3 


45 

. 6 

62 

. 6 








Figure 3. (concluded) Output for sample case. 
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Section 4 

Source Code Listing 


The remainder of this volume gives the listing of the routines shown in figure 1. The 
routines can be categorized as follows. The first group compises new code as described in 
Volume 1. Then there are two routines taken verbatum from the Smith code: DSWK and 
WAVE. Finally, there is a series of routines for inverting the (double precision, real) matrix of 
influence coefficients. These were taken from LINPAC (see LINPAC User’s Guide, SIAM, 
Philadelphia, 1979) and are not shown here since they are in the public domain and are 
commonly available. Of course, they are included on the disk with the rest of the source code. 
The routine MATINV is included, since this was written for the present purposes to call the 
LINPAC routines. 



o o 


PROGRAM CUP 2D 


c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 
c . . 

c . . 


Calculates the unsteady loading and associated acoustic and vorticity waves 
on 2 mutually interacting blade rows in a 2 dimensional, linear, subsonic 
analysis. Blade rows can be rotor/egv, igv/rotor, or rotor/rotor, depending 
input values of rotor rotational Mach numbers, MY1 and MY2. 

Simultaneous solution for flow tangenoy for 2 blade rows, NH harmonics, and 
NP panels via inversion of the matrix coupling equation KMATRX * LOAD- WAS H . 

A disturbance upwash distribution at either, or both, blade rows is generated 
either from direct user input or from Silverst.ein ' s wake formulas. The code 
code then finds the unsteady loading (LOAD) that produces an upwash (WASH) 
that just cancels the disturbance wash. These loads are used to find the 
acoustic waves and the sound power. Blade row self-effect is computed via 
subroutines from Smith Code. Effect of each row on the other is computed 
via an extension of Smith's theory by D.B. Hanson. Reflections at inlet and 
exit interfaces are treated with reflection and transmission coefficients 
derived from continuity of mass & momentum based on an actuator disk model. 
Overall theory documented in NASA CR-4506 , Volume I. Comments in this listing 
refer to equation numbers in the same Contractor Report. Coding by Hanson. 

This routine is the main program. 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 


DOUBLE PRECISION MX A , MX B , MXC , MY 1 , MY 2 , MS , LAM 1 , LAM2 , 

> KMATRX (1020, 1020) 

DOUBLE COMPLEX WREXT ( 5, f » 1 ) , WSEXT ( 5, 51 ) , LR(5,51), LS(5,51), 

> ALF (9, -5:5, -5:5) , 


> KRUP ( 5 , -5 

> R 1 2 (-5:5, 

> T14 (-5:5, 
INTEGER Bl, 


5,51) ,KRDN(5, -5:5,51) ,KSUP(5, -5:5,51) , KSDN ( 5 , -5 : 5 , 51 ) , 
5:5), R21 (-5:5, -5:5) , R1 3 ( -5 ; 5 , - 5 : 5 ) , R3 1 { -5 : 5 , -5 : 5 ) , 
5:5) , T2 8 ( - 5 : 5 , -5:5) , T 3 8 ( - 5 : 5 , -5:5) 

B2 , BETA ( -5 : 5 , - 5 : 5 ) 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
WRITE ( * , * ) ' ' 

WRITE {*,*) ' Code CUP2D for coupled cascade aeroacoustics - Version 

> 1.0 ' 

WRITE (*,*) 'Developed for NASA-Lewis by Pratt & Whitney under Cont. 
>ract NASH-25952 - Task 10' 

WRITE (*,*)' Theory documented in NASA CR-4506, Volume I' 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c.. Read input, generate wavenumbers and reflection coefficients, establish 
c.. external disturbance vectors from wake formulas or direct input, and 
c.. print input geometry and flow conditions. 

CALL READIN { NH , NP, Bl , B2 , Cl , C2 , SCI , SC 2 , CT1 , CT2 , ST1 , ST2 , XS , XA, XC, 

> MXA, MXB, MXC, MS , MY1 , MY2 , LAM1 , LAM2 , ALF, BETA, R12 , R2 1 , R1 3 , R3 1 , 

> T14,T28,T38, AA, AB, AC, WREXT, WSEXT , POPSA, POPSC ) 

c.. Generate matrices of influence functions 

CALL INFFNS ( NH , NP , Bl , B2 , SCI , SC2 , CT1 , CT2 , ST1 , ST2 , XS , 

> MXB, LAM1, LAM2 , ALF, BETA, R12, R21 , R13, R31 , T14 , T28, T38, 

> KMATRX, KRUP , KRDN, KSUP , KSDN ) 

c.. Solve coupled system of equations for loading by matrix inversion. 

CALL SOLVE ( NH , NP , KMATRX , WREXT , WSEXT , LR , LS } 

c.. Compute output waves and print their sound pressure and sound power. 

CALL OUTPUT (NH , NP, Bl , B2 , MXA , MXB , MXC, MS, MY1 , MY2 , ALF, AB/AA, AB/AC, 

> XS, XA, XC, POPSA, POPSC, KRUP, KRDN, KSUP, KSDN, LR , LS ) 


END 


Figure 4. Source code for CUP2D. 
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c . . 
c . . 
c . . 
c . . 


c . . 


c . . 


c . . 
c . . 
c . . 
c . . 
c . . 


c . . 


c . . 


c . . 


SUBROUTINE READIN <NH , NP , B1 , B2 , Cl , C2 , SCI , SC2 # CT1 , CT2 , ST1 , ST2 , 

> XS, XA, XC,MXA,MXB,MXC,MS,MY1,MY2, LAM1,LAM2, ALF, BETA, 

>R12 , R21 , R13 , R31 , T14 , T28 , T38 , AA, AB, AC, WREXT, WSEXT, POPSA, POPSC) 

Reads input from data set on disk {UNIT 8); generates constants for the 
Smith common block; and calls routines for alpha and beta wavenumbers, 
reflection and transmission coefficients, and external input velocity 
vectors . 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION MXA , MXB , MXC , MS , MY 1 , MY2 , LAM1 , LAM2 
INTEGER Bl, B2, BETA ( - 5 : 5 , - 5 : 5 ) 

DOUBLE COMPLEX WREXT ( 5 , 51 ) , WSEXT ( 5 , 51 ) , ALF ( 9 , -5 : 5 , - 5 : 5 ) , 

R12 (-5:5, -5:5) ,R21 ( -5 : r > , -5 ; 5) ,R13 (-5:5, -5:5) ,R31 (-5: 5, -5: 5) , 
T14 (-5:5, >5:5) ,T28 ( -5 : 5, -5: 5) ,T38 (-5: 5, -5 :5) 

CHARACTER*? 0 COMMENT 


Read and compute data for common block 
OPEN (UNIT=8, FI LE= ' cupin.dat ' ) 

Read comment 

READ (8,*) COMMENT 

Read geometry from disk file and compute normalized chords 
READ { 8, * ) Bl, B2 , SCI, SC2 , XS 

Cl = 6 . 28318S3D0/ ( B1*SC1 ) ! Chord/radius , front row 

C2 = 6 . 2831 853D0 / { B2*SC2 ) ! Chord/radius, back row 


Read number of panels and number of harmonics 
READ ( 8 , * ) NP, NH 

Read operating conditions from disk file. 1 & 2 refer to upstream and 
downstream blade rows. A, B, & C refer to regions upstream of blade row 1, 
between blade rows, and downstream of blade row 2. POPSA & POPSC = ambient 
pressure in upstream and downstream regions divided by sea level standard 
pressure . 

READ ( 8 , * ) MXA, MXB, MXC, MS ! Axial & swirl Mach # # s 

READ (8,*) AA, AB, AC, RHOA, RHOB, RHOC ! Sound speeds & densities 

READ (8,*) MY1 , MY 2 ! Blade rotational Mach #'s 

Read x locations for pressure output, (x over radius from front row LE.) 

READ (8,*) XA, XC ! Make XA < 0 and XB > 0 . 


Compute remaininq 
ST1 = ( -MY1 +MS) 
ST2 = ( -MY2+MS) 
CT1 = SQRT { IDO 
CT2 = SQRT {IDO 


items for passaqe to other 
/ SQRT ( MXB* * 2 +~ { -MY1+MS) * * 
/ SQRT ( MXB* *2 4 { -MY2+MS) * * 
ST1 * *2 ) 

ST2* *2 ) 


routines 

2) ! Sine(thetal) 

2) ! Sine(theta2) 

! Cosine { thetal } 
! Cosine ( theta2 ) 


LAM1 = B2* (MY1-MY2) / (MXB/CT1) 
LAM 2 = Bl* (MY1-MY2) / (MXB/CT2) 


! Reduced freq, front row 
! Reduced freq, back row 


Compute axial wavenumbers (alpha's) and tanqential wavenumbers (beta's) 
CALL ALFBET { Bl , B2 , NH , MXA , MX B , MXC , MS , MY1 , MY2 , 

> AA, AB, AC, ALF, BETA) 


Get reflection and transmission coefficients 

CALL RTCOEF { NH , Bl , B2 , ALF , BETA , AA , AB , AC , RHOA , RHOB , RHOC , 

> MXA , MXB , MXC , MS , MY 1 , MY 2 , R12 , R21 , R13 , R31 , T14,T2 8,T3 8) 

Read wake input data and compute upwash vectors. Leave unit 8 open to 
read from GTWAKE. 

CALL GTWAKE ( NH , NP , Cl , C2 , SCI , CT1 , CT2 , ST1 , ST2 , XS , INTY PE , CD, 

> WREXT, WSEXT) 


CLOSE (8) 


Figure 4. (continued) Source code for CUP2D. 
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c.. Compute Pambient./Pstandard to be used in OUTPUT for SPL's 
POPSA = RH0A*AA**2/ (1 .4*32 .2*2116) 

POPSC = RH0C*AC**2/{1. 4*32. 2*2116) 

c,. Print input data 

CALL PRNTIN ( NH , NP , B1 , B2 , Cl , C2 , SCI , SC2 , CT1 , CT2 , ST1 , ST2 , XS , 

> AA , AB, AC , RHOA , RHOB , RHOC , MXA , MXB , MXC , MS , MY1 , MY 2 , 

> LAM1 , LAM 2 , COMMENT, INTYPE, CD, WREXT, WSEXT, POPSA, POPSC ) 

RETURN 

END 


Figure 4. (continued) Source code for CUP2D. 
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SUBROUTINE ALFBET (B1 , B2 , NH, MXA, MXB, MXC, MS, MY1 , MY2 , 

> AA, AB, AC, ALF, BETA) 

c.. Computes alpha and beta wavenumbers from formulas derived in appendix B. 
c.. The alphas are Smith's normalized by source radius R rather than by chord, 
c.. Prints message if a resonance condition occurs for any combination of n & k. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION MXA , MX B , MXC , MS , MY 1 , MY2 
INTEGER Bl,B2,BETA(-5:5,-5:S) 

DOUBLE COMPLEX ALF {9, -5:5, -5:5) 

DENOMA = l.ODO - MXA* *2 

DENOMB = l.ODO - MXB* *2 

DENOMC = l.ODO - MXC* *2 

DO 100 N = -NH, NH 
DO 100 K = -NH, NH 

IF ( ABS (N ) + ABS { K ) . EQ . 0) GOTO 100 

BETA ( N, K ) = - (N*B1 - K*B2) 

OMEGA = N*B1*MY1 - K*B2 *MY 

EA = DENOMA* BETA (N, K) **2 - ( OMEGA *AB/AA } **2 

EB = DENOMB* BETA (N, K) * *2 - (OMEGA + BETA (N, K) *MS) * *2 

EC = DENOMC* BETA (N, K) **2 - ( OMEGA *AB/ AC ) * *2 

FA = MXA*OMEGA*AB/AA 

FB = MXB* (OMEGA + BETA ( N, K ) *MS) 

FC = MXC* OMEGA *AB /AC 

c.. Check for resonance in upstream region, swirl region, and downstream region 
c.. E=0 for resonance, E < 0 for propagation, E > 0 for decay. Then compute 
c.. alphas for pressure waves. Use Eqs . B-8,9,10 in Region B and variations 
c.. for Regions A and C. 

c. . Do Region A first (upstream of swirl region) . 

IF (EA .EQ. 0.0) THEN 
WRITE ( * , 1 ) N , K 

1 FORMAT ( IX, ' Resonance in Region A for N=',I3,' K=',I3) 

STOP 'Execution terminated due to resonance' 

ELSE IF (EA .LT. 0.0) THEN 

ALF ( 4 , N , K ) = ( FA + SIGN ( SQRT { - EA ) , FA) ) /DENOMA 

ALF { f>, N, K) = (FA - SIGN { SQRT ( -EA ) , FA) ) /DENOMA 

ELSE 

ALF ( 4 , N, K) = CMPLX (FA, -SQRT ( EA) ) /DENOMA 
AL F ( T> , N , K ) = CMPLX ( FA, + SQRT ( EA) ) /DENOMA 
ENDIF 

c.. Do Region B next (swirl region). 

IF (EB .EQ. 0.0) THEN 
WRITE ( * , 2 ) N , K 

2 FORMAT ( IX, ' Resonance in Region B for N=',I3,' K=',I3) 

STOP 'Execution terminated due to resonance' 

ELSE IF (EB .LT. 0.0) THEN 

ALF ( 1 , N, K ) = (FB + SIGN ( SQRT ( -EB ) , FB) ) /DENOMB 

ALF ( 2 , N, K ) = (FB - SIGN ( SORT ( -EB ) , FB) ) /DENOMB 

ELSE 

ALF ( 1 , N, K } = CMPLX (FB, -SQRT ( EB )) /DENOMB 
AL F ( 2 , N , K ) = CMPLX (FB, + SQRT ( EB) ) /DENOMB 
ENDIF 

Finally, do Region C (downstream of swirl region) 

IF (EC .’EQ. 0.0) THEN 
WRITE ( * , 3 ) N , K 

FORMAT ( IX, ' Resonance in Region C for N=',I3,' K=',I3) 

STOP 'Execution terminated due to resonance' 

ELSE IF (EC .LT. 0.0) THEN 

ALF ( 7 , N , K ) = (FC + SIGN ( SQRT ( -EC) , FC) ) /DENOMC 


Figure 4. (continued) Source code for CUP2D. 
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ALF ( 8 , N, K) = (FC - SIGN ( SQRT ( -EC) , FC) ) /DENOMC 
ELS E 

"aLF(7,N,K) = CMPLX(FC, -SQRT ( EC) } /DENOMC 
ALF ( 8 , N, K) = CMPLX(FC, + SQRT ( EC) ) /DENOMC 
ENDIF 

c.. Wavenumbers for vorticity waves in Region B from eq. B-ll and variations 
c.. for Regions A and C. 

ALF ( 3 , N, K ) = -(OMEGA + MS * BETA ( N, K ) ) /MXB 
ALF ( 6 , N, K ) = - OMEGA /MXA 
ALF ( 9 , N, K } = - OMEGA /MXC 

100 CONTINUE 

RETURN 

END 


Figure 4. (continued) Source code for CUP2D. 
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SUBROUTINE RTCOEF(NH, B1 , B2 , ALF, BETA , AA, AB , AC, RHOA , RHOB , RHOC, 

> MXA , MXB , MXC , MS , MY I , MY 2 , R12 , R21 , R13 , R31 , T14,T28,T38) 

c-- This subroutine computes reflection and transmission coefficients for the 
c.. inlet and exit for the transverse velocity component. Coefficients derived 
c.. in NAS’A CR-4506, Volume I, appendix D. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION MXA, MXB, MXC, MS , MY1 , MY2 
INTEGER B1 , B2 , BET , BETA ( -5 : 5 , - 5 : 5 ) 

DOUBLE COMPLEX Cl , C2 , C3 , C4 , C8 , C9 , FI , F2 , F3 , F4 , F8 , F9 , 

> G1 , G2 , G3 , G4 , G8 , G9 , EO, El, E2 , E3 , E4 , E5, E6, E7 , E8, E9, E10, 

> Ell , E12 , E13 , E14 , E15 , E16 , E17 , E18 , E19 , E20 , E21 
DOUBLE COMPLEX ALF( 9, -5 : 5, -5 : 5) , 

>8121-5:5,-5:5), R21 { -5 : 5, -5 : 5) , R13 { -5 : 5, -5 : 5) , R31 (-5:5, -5: 5) , 

> T14(-5:5,-5:5) , T28 ( -5 : 5 , -5 : 5) , T38 ( -5 : 5, -5 : 5) 

c.. Compute rho*cO for 3 reqions 
ROCA = RHOA*AA 
ROCB = RHOB* AB 
ROCC = RHOC* AC 


C. . Compute reflection and transmission coefficients 
DO 10 N = -NH, NH 
DO 10 K = -NH , NH 

IF ( (ABS(N)+ABS(K) ) . EQ . 0) GO TO 10 

BET = BETA (N, K) 

OMEGA = N* B1 *MY1 - K*B2*MY2 


C.. Coefficients from continuity equations 
Cl = ROCB* ( ( 1D0-MXB**2 ) *ALF( 1 , N, K) - 
C2 = ROCB* ( ( 1 DO- MXB* *2 ) *ALF(2, N, K) - 
C3 = ROCB* {-BET) 

C4 = ROCA* ( ( 1D0 -MXA* *2}*ALF(4,N, K) - 
C8 = ROCC* ( { IDO -MXC* *2 ) * ALF ( 8 , N, K ) - 
C9 = ROCC* {-BET) 


MXB* (OMEGA+MS*BET) ) 
MXB* ( OM EGA+ MS* BET ) ) 

MXA* OMEGA* AB/AA ) 
MXC*OMEGA*AB/AC ) 


C.. Coefficients from axial momentum equations 

FI =ROCB* { { 1D0+MXB* *2 ) * ( OMEGA + MS* BET) - ( 1D0-MXB**2) *MXB*ALF { 1 , N, K ) ) 
F2=ROCB* { { 1D0 + MXB* *2 ) * { OMEGA+MS* BET) - ( 1D0-MXB* *2 ) *MXB*ALF { 2 , N, K ) ) 
F3=ROCB* ( 2D0*MXB*BET) 

F4=RC>CA* ( { 1D0+MXA* *2 ) * OMEGA - ( 1D0 -MXA* *2 ) * MXA* ALF { 4 , N, K) *AA/AB } 
F8 = ROCC* ( ( 1 DO + MXC* * 2 ) * OMEGA - ( 1D0 -MXC* *2 ) * MXC* ALF ( 8, N, K) *AC/AB) 
F9=ROCC* { 2D0 *MXC*BET*AC/AB ) 


C.. Coefficients from transverse momentum equations 

Gl = ROCB* ( MXB* BET- MXB* MS* ( OMEGA+MS * BET ) + (1D0-MXB**2) *MS*ALF(1,N, K) ) 
G2 = ROCB* ( MXB* BET- MXB* MS* (OMEGA+MS* BET) + (1D0-MXB**2) *MS*ALF(2,N, K) ) 
G3= ROCB* (MXB* ALF (3 , N, K ) -MS* BET) 

G4 = ROCA* (MXA* BET* AA/AB) 

G8= ROCC* (MXC* BET* AC/AB) 

G9= ROCC* (MXC* ALF (9, N, K) *AC/AB) 


EO 

= 

C4 * FI 

- Cl * F4 

El 

- 

C4*F2 

- C2*F4 

E2 

= 

C3*F4 

- C4*F3 

E3 

= 

C4*G1 

- Cl *G4 

E4 

- 

C4*G2 

- C2*G4 

E5 

= 

C3*G4 

- C4*G3 

E6 

= 

C2*G1 

- Cl *G2 

E7 


C3*G2 

- r*2*(;3 

EH 

— 

C2*F1 

- Cl * F2 

E9 


C 3 * F 2 

- C2*F3 

E10 

= 

C1*FH 

- c:b*fi 


Figure 4. (continued) Source code for CUP2D. 
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E12 = C1*G8 - C8*G1 
El 3= C9*G1 - Cl *G9 
E14 = C8*F2 - C2*F8 
El 5= C8*F9 - C9*F8 
El 6= C8*G2 - C2 *G8 
E17= C8*G9 - C9*G8 
El 8= C8*F3 - C3*F8 
El 9= C8*G3 - C3*G8 
E2 0= C1*F3 - C3*F1 
E21= C1*G3 - C3 *G1 


C.. Reflection co 
R1 2 { N , K ) = 
R 1 3 { N , K ) = 
R21(N,K) = 
R 3 1 { N , K ) = 

C. . Transmission 
T14(N,K) = 
T2 8 { N, K ) = 
T3 8 ( N, K } = 

10 CONTINUE 


efficients 
( E2 *E3- E0*Ef>) / 

( El *E3- E0*E4 ) / 

(E15*E16-E14*E17) / 
(E15*E19-E17*E18) / 

::oef f icients 
(E7*E8-E6*E9) / 

(E6*E11-E8*E13) / 

{ E11*E21-E13*E20) / 


(El*E r >-E2*E4) 

( El *ES- E2 *E4 ) 

(E12*E15-E10*E17) 

(E12*E15-E10*E17) 


( E4 * E9- El *E7 ) 

( El 0* El 3 -Ell* El 2) 
(E11*E12-E10*E13) 


RETURN 

END 

C 

C 


Figure 4. (continued) Source code for CUP2D. 



SUBROUTINE GTWAKE (NH, NP, Cl , C2 # SCI , CT1 , CT2 , ST1 , ST2 # XS # 

> INTYRE, CD, WREXT, WSEXT) 

c.. This routine reads data from the input file and generates the upwash vectors 
c.. WREXT and WSEXT, representing external excitation of the system. The 
c.. disturbance can be described via various optional methods specified by the 
c.. input value of INTYPE as follows. 

c.. INTYPE = 1 is used to represent viscous wakes via the Silverstein formulas, 
c.. The wake is specified by the drag coefficient CD, i.e. by only one number, 
c.. The upwash vectors are then computed from the formulas in appendix E. 
c.. INTYPE = 2 is the same as INTYPE 1 except that the wake harmonics above BPF 
c.. are set to zero. 

c.. INTYPE = 3 is the same as INTYPE 1 except, that the velocity defect harmonics 
c.. are input by the user rather than being computed from the wake formulas, 
c.. For INTYPE = 4 the real and imaginary parts of WREXT and WSEXT are simply 
c.. read from the file. Here the upwash vectors are completely specifed by the 
c . . user for all harmonics: N = 1...NH and all control points I = 1...NP. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION FW(5) 

DOUBLE COMPLEX EXPON, WREXT ( 5 , 51 } , WSEXT ( S , r >l ) , AI 

PI = 3.141 S926SD0 

AI = (0. 0D0, 1 .01)0) 

c.. Read INTYPE to identify type of input disturbance 
READ ( 6 , * ) INTYPE 
IF (INTYPE . EC) . 1) THEN 

c.. For INTYPE=1, read drag coefficient and use Silverstein formulas as given 
c.. in appendix E. 

READ (8,*) CD ! Draq coefficient 

DO 34 N = 1, NH 
DO 34 I = 1, NP 

221 = 0 . f>D0* ( IDO -COS ( PI * ( 2 . 0D0* I - 1D0 ) / ( 2 . 0D0*NP ) ) ) 

Z1 = (XS/C1 4 Z2I/ (C1/C2 ) *CT2) /CT1 
WCOW1 = 1 .21D0*SQRT(CD) / (Z1 - 0.70D0) 

Y0C1 “ 0. 68D0*SQRT(CD* (Z1 - 0.8f>D0)) 

0 = 1 ,7724r.DO*SCl/YOCl*CTl 

FN = 3 . r >4 4 91D0/Q* WCOW1 * EXP ( - ( PI *N/0) * * 2 ) 

EXPON =EXP (AI *2D0*PI *N* ( 

> (CT2*ST1/CT1-ST2)MC2/C1)/SC1*Z2I + (XS/C1 ) /SCI *ST1/CT1 ) ) 

WSEXT ( N , I ) = CT2/CT1* (ST2*CT1 -CT2*ST1 ) *FN*EXP0N 

WREXT {N , I ) = (0D0, 0D0) 

34 CONTINUE 

c.. For INTYPE=2, read drag coefficient and use Silverstein formulas as given 
c.. in appendix E. BUT, for harmonic order > BPF, set upwash to zero. 

ELSE IF (INTYPE . EQ . 2) THEN 

READ(8,*) CD ! Draq coefficient 

DO 24 N = 1, NH 
DO 24 I = 1, NP 

Z2I = 0 . f>D0* ( IDO -COS (PI * ( 2 . ODO* I - IDO ) / ( 2 . 0D0*NP) ) ) 

Z1 = (XS/C1 4 Z2I/ (C1/C2) *CT2 } /CT1 
WCOW1 = 1 .21D0*SQRT(CD) / (Z1 - 0.7 ODO) 

YOC1 = 0.68D0*SQRT(CD*(Z1 - O.ftSDO}) 

Q = 1 .77245DO*SC1/YOC1*CT1 

FN = 3.f»4491DO/Q*WCOWl*EXP(- (PI*N/Q) **2) 

EXPON =EXP(AI # 2D0^PI*N* ( 

> (CT2 *ST1 /CT1 -ST2 ) * (C2/C1 ) /SCI *Z2 1 4 (XS/C1 ) /SCI *ST1/CT1 ) ) 

WSEXT (N, I) = CT2/CT1* ( ST2 *CT1 -CT2 *ST1 ) *FN*EXPON 

IF (N .G T. 1) WSEXT (N, I) = (ODO, ODO) 

WREXT (N, I) = (ODO, ODO) 

24 CONTINUE 


Figure 4. (continued) Source code for CUP2D. 
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c.. For INTYPE=3, read harmonics of an upwash that convects with the mean flow, 
c,, This is just like INTYPE 1 above except that FW(N) is read from input here 
c. . and is independent of x. By contrast, for INTYPE = 1 above, FN is computed 
c.. from the Silverstein formulas as a function of chordwise position on the 
c.. downstream blade row. To interpret these formulas, see appendix E. 

ELSE IF (INTYPE . EQ. 3) THEN 

READ ( 8 , * ) (FW(N) ,N=1,NH) 

DO 14 N = 1 , NH 
DO 14 I = 1, NP 

Z2I = 0 . f>D0* ( IDO - C0S1PI* (2. 0D0M-1D0) / (2 .0D0*NP) ) ) 

EXPON =EXP(AI*2D0*PI*N* ( 

(OT2 *ST1 /CT1 -ST2 ) * (C2/C1 ) /SC1*Z2I + (XS/C1 ) /SCI *ST1 /CT1 ) ) 

WSEXT (N, I ) = CT2 /CT1 * ( ST2*CT1 -CT2 *ST1 ) *FW (N) * EX PON 
WREXT ( N , I ) = (0D0, 0D0) 

14 CONTINUE 

c.. For INTYPE=4 , read real and imaginary parts of vectors of external velocity 
c.. disturbance as direct input. Rotor input first, then stator. 

ELSE IF (INTYPE . EQ . 4) THEN 
DO 10 N=1 , NH 
DO 10 1=1, NP 

READ (8,*) WREXTR, WREXTI 
WREXT (N, I) =CMPLX (WREXTR , WREXTI ) 

10 CONTINUE 

DO 12 N=1 , NH 
DO 12 1=1, NP 

READ (8,*) WSEXTR,WSEXTI 
WSEXT (N, I ) =CMPLX (WSEXTR, WSEXTI ) 

12 CONTINUE 


ELSE 

STOP 'Input, type (INTYPE) for WSEXT AND WREXT not defined' 
END IF 

RETURN 

END 


Figure 4. (continued) 
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SUBROUTINE PRNTIN (NH, NP , B1 , B2 , Cl , C2 , SCI , SC2 , CT1 , CT2 , ST1 , ST2 , XS, 

> AA , AB , AC , R HOA , RHOB , RHOC , MXA , MXB, MXC, MS , MY 1 , MY2 , 

> LAM1,LAM2, COMMENT, INTYPE, CD , WREXT , WSEXT, POPSA, POPSC) 

c.. This routine prints the input data (some of it manipulated) to the screen. 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION MXA , MXB , MXC , MS , MY1 , MY2 , LAM1 , LAM2 
INTEGER Bl, B2 

DOUBLE COMPLEX WREXT(5,51), WSEXT (5, 51) 

CHARACTER *70 COMMENT 

c.. Print input data and computed quantities 
WRITE (*,*) ' ' 

WRITE)*,*) 'COMMENT: COMMENT 

CALL TIMDAT ! WRITES TIME AND DATE OF EXECUTION 

WRITE!*,*) ' ’ 

WRITE (*,11) Bl, B2 

11 FORMAT ( IX , ' Bl= ' , 13 , ' B2=',I3) 

WRITE!*, 13) SCI, SC 2 

13 FORMAT (IX, 'Gap/Chord (1) = ’,F6.3,', Gap/Chord (2) = ',F6.3) 

WRITE (*,31) XS 

31 FORMAT ( IX, ' (Rotor LE to Stator LE)/ (Local Radius of Rotor ) = ’ , F6 . 3 , 

>' ( input) ' ) 

WRITE!*, 35) ( XS-C1 *CT1 ) /Cl 

35 FORMAT ( IX ,' Axi al Spacing Between Blade Rows/Rotor Chord ' , F7 . 4 , ' 

> (computed) ') 

IF (INTYPE .LT. 3) WRITE!*, 1) CD 
1 FORMAT (IX, 'Drag Coefficient = ',F5.3) 

WRITE!*,*) ' ' 

WRITE!*, 15) NP, NH 

15 FORMAT! IX, 'Number of panels= ' , 12 , ' , Number of harmonics= ' , 12 ) 

WRITE!*,*) ' ' 

WRITE!*,*)' Mxa Mxb Mxc Ms Myl My2 Mrell Mre 

>12 ' 

WRITE!*, 17) MXA, MXB ,MXC, MS, MY1 , MY2 , MXB/CT1, MXB/CT2 
17 FORMAT ( IX , 8F7 . 3 ) 

WRITE!*,*) ' ' 

WRITE!*,*) ' RHOa RHOb RHOc Aa Ab Ac ' 

WRITE!*, 19) RHOA, RHOB, RHOC, AA, AB, AC 

19 FORMAT ( IX , 3F8.5, 3F8.1) 

WRITE!*,*) ' ' 

WRITE!*,*) 'Remainder of printout is computed from input above' 

WRITE!*, 21) LAM1*C1 , LAM2*C2 

21 FORMAT (IX, 'Smiths reduced freqs @ BPF front row, rear row =',2F8.3) 
WRITE!*, 27) 57 . 2957R*ASIN ( ST1 ) , 57 . 2957 8*ASIN (ST2 ) 

27 FORMAT (IX, 'Theta 1 , Theta2 (in degrees) = ', 2F8.3) 

WRITE ( * , 29 ) 57 . 29578*ATAN (MS/MXB) 

29 FORMAT! IX, ' Swirl Angle (in degrees) = ', F7.2) 

WRITE!*, 16) POPSA, POPSC 

16 FORMAT (IX, 'Ambient Pressure /Sea Level STD (uptream, downstream) =' 

> , 2F7 . 3 ) 

WRITE!*, *) ' ' 

WRITE!*,*) ' EXTERNAL VELOCITY IMPOSED ON CASCADE’ 

WRITE!*,*) ' N I WREXT ( real , imag) WSEXT(real, imag)' 

DO 20 N=1 , NH 
DO 22 1=1, NP 

WRITE!* , 33) N, I , WREXT (N, I ) , WSEXT (N, I) 

33 FORMAT (IX, 2 14, ' \2FR.4,' ',2F8.4) 

22 CONTINUE 

WRITE!* , * ) ' ' 

20 CONTINUE 

WRITE!*,*) ' ’ 

RETURN 

END 

C 


Figure 4. (continued) Source code for CUP2D. 
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SUBROUTINE TIMDAT ! This routine is specific to the Sun Workstations 

CHARACTER*24 FDATE ! Modify this routine for other computers 

WRITE (*,*) 'Time of execution: \ FDATE (} 

END 


Figure 4. (continued) Source code for CUP2D. 
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SUBROUTINE INFFNS (NH, NP, B1 # B2 , SCI , SC2, CT1 , CT2 , ST1 , ST2 , XS, 

> MXB, LAM , LAM2 , ALF , BETA, R12, R21, R13, R31 , T14 , T28, T38, 

> KMATRX , KRUP, KRDN, KSUP , KSDN ) 

c.. Calls routines to compute elements of KMATRX, the matrix of influence 
c.. functions. KMATRX is then returned to the main program. Algebra is 
c.. based on NASA CR-4506, Volume I, Section 3.3. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION MXB , LAM1 , LAM 2 , KMATRX ( 1 02 0 , 1 020 ) 

INTEGER Bl, B2 , BETA ( -5 : 5, -5 : 5 ) 

DOUBLE COMPLEX ALF (9, -5:5, -5 : 5) , 

> R 1 2 { - 5 : 5 , - 5 : 5 ) , R21 { -5 : 5 , -5 : 5 ) , R13 ( -5 : 5 , -5 : 5) , R31 ( -5 : 5, -5: 5) , 

> T14{-5:5,-5:5) , T28 ( -5 : 5, -5 : 5 ) , T38 ( -5:5, -5: 5) , 

KRUP {5,-5 : 5, 51} , KRDN ( 5 , - 5 : 5 , 51 ) , KSUP ( 5 , - 5 : 5 , 51 ) , KSDN (5, -5: 5, 51 ) 

Cl = 6 . 2831 853D0/ (B1*SC1 ) ! chord/ radius , row # 1 

C2 = 6 . 2831 85300/ ( B2*SC2 ) ! chord/ radi. us , row # 2 


C.. Zero the Kmat.rix before starting to compute the elements 
DO 10 MU = 1, 4 *NP*NH 
DO 10 NU = 1, 4 *NP*NH 

KMATRX {MU, Nil) = 0 . 0D0 
10 CONTINUE 

C.. Effect of rotor loading on stator 

CALL RCOEFS ( NH , NP , Cl", C2 , SCI , CT1 , CT2 , ST1 , ST2 , XS , 

> MXB, LAM1 , ALF, BETA, R12 , R21 , R1 3 , R3 1 , T1 4 , T2 8 , T3 8 , KMATRX , KRUP, KRDN) 

C.. Effect of stator loading on rotor 

CALL SCOEFS ( NH , NP , Cl , C2 , SC2 , CT1 , CT2 , ST1 , ST2 , XS , 

> MXB, LAM2 , ALF, BETA, R12 , R2 1 , R1 3 , R31 , T14, T2 8 , T3 8 , KMATRX , KSUP, KSDN) 

C.. Effect of rotor loading on rotor 

CALL GENKRR {NH, NP, Cl] C2 , SCI , SC2 , CT1 , ST1 , MXB , LAM1 , KMATRX ) 

C,. Effect or stator loading on stator 

CALL GENKSS (NH , NFL Cl , C2 , SCI , SC 2 , CT2 , ST2 , MXB , LAM2 , KMATRX ) 

RETURN 

END 


Figure 4. (continued) 
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SUBROUTINE RCOEFS (NH, NP, Cl , C2 , SCI ,CT1 , CT2 , ST1 , ST2 , XS, 

> MXB, LAM1 , ALF , BETA , R1 2 , R21 , R1 3 , R3 1 , T14 , T28 , T38 , KMATRX, KRUP, KRDN ) 

c.. Generates elements of the matrix of influence coefficients KMATRX that give 
c.. the upwash caused by rotor loading at control points on the stator and 
c.. rotor. These are computed from KRS (N, K, I , J) , effect of rotor on stator, 
c.. and KRR ' (N, I , J ) 1 , effect of stator on stator. , The prime on KRR ' indicates 
c.. that only the part of KRR associated with waves reflected from the actuator 
c.. disk is computed here. The remainder is computed in GENKRR using original 
c.. routines from the Smith code. N counts the rotor loading harmonics, I the 
c.. control points, and J the load elements. K counts the cascade wave index in 
c. . the rotor frame which becomes the time harmonic index in the stator frame. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION MXB, MR1 , MR2 , LAM1 , KMATRX (1020, 1020) 

INTEGER BET, BETA (-5:5, -5:5) 

DOUBLE COMPLEX AI , ALF1 , ALF2 , ALF3 , EXPE21 , EXPE3 1 , EXPE1 , 

> El I , E2I , E3 I , KR1,KR2, KR3, V1,V2,VR1, VR2, VR3 , KRS, KRR, 

> ALF(9, -5:5, -5:5) , KRUP ( 5 , - 5 : 5 , 51 ) , KRDN ( 5 , *5 : 5 , 51 ) , 

> R12(-5:5,-5:5) , R21 ( -5 : 5, -5 : 5) , R13 ( -5 : 5, -5: 5) , R31 (-5 :5, -5:5) , 

> T14(-5:5,-5:5) , T28 ( -5 : 5, - 5 : 5) , T38 ( -5 : 5 , -5 : 5) 

AI = (0.0D0, 1.0D0) 

PI = 3.141 59265D0 

WRITE ( * , * ) 'Entering RCOEFS ' 

MR1 = MXB/CT1 

MR2 = MXB/OT2 

CON = MR1/SC1 

XE = XS 4 C2*C T2 

NPNH2 = NP*NH*2 
DO 10 N = 1, NH 

DO 10 K = -NH, NH 
BET = BETA { K , N } 

ALF1 = ALF ( 1 , K , N) 

ALF2 = ALF ( 2 , K, N) 

ALF3 = ALF ( 3 , K, N) 

EXPE21 = EXP (AI* (ALF2-ALF1 ) *XE) 

EXPE31 = EXP ( AI * ( ALF 3- ALF 1 ) *XE } 

EXPE1 = EXP ( -AI *ALF1*XE) 

c.. Get Smith's factors for transverse velocity components 
CALL GETVS ( BET , N* LAM 1 ,MXB, DTI , CT1 , VI , V2 , V3 } 

DO 10 J = 1, NP ! Loop on load elements 

Z0J = 0 . 5D0M 1D0 * COD ( PI* (.7- 1D0} /NP) ) ! Chordwise locations for loads 

KR 1 = CON*V1*EXP(-AIMALF1*CT1+BET*ST1)*C1*ZOJ) ! 

KR2 = CON*V2*EXP(-AI*(ALF2*CT1+BET*ST1)*C1*ZOJ) ! eq . 43 

KR3 = CON*V3*EXP(-AIMALF3*CT1 + BET*ST1)*C1*ZOJ) ! 

c.. Compute VR1, VR2, VR3 from eq . 44. 

c.. Hold EXP ( - AI *ALF1 *XE) out of VR 1 now to avoid overflow in KRS later 
VR1 = ( (KR1*R12 (K , N) +KR2) *R21 (K,N) *EXP (AI *ALF2 *XE) 

> 4 ( K R 1 * R 1 3 ( K , N ) + KR 3 ) * R 3 1 (K,N) * EXP ( AI *ALF3*XE) ) 

> / (1D0-R12 (K,N)*R21 ( K , N} *EXPE21-R1 3 ( K , N) *R31 (K,N) *EXPE31 ) 

VR2 = (KR1+VR1*EXPE1) *R12(K,N) 

VR3 = { KR 1 4 VR 1 * EXPE1 ) *Ri 3 ( K , N } 

c.. KRUP and KRDN are passed out of the subroutine for later use in computing 
c.. pressure in the non-swirl regions upstream and downstream. 

KRUP { N, K , J ) = (KRI 4 VR 1 * EX PEI ) * T14(K,N) ! eq . 85 

KRDN { N , K , J ) = (KR 24 VR 2 ) *T28 (K,N) *EXP ( AI *ALF2*XE) ! eq. 97 

> 4 (KR3+VR3) *T3B (K,N) *EXP(AI *ALF3*XE) 


Relative Mach #, row 1 
Relative Mach #, row 2 
A constant 

Axial coordinate of stator exit 
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c.. Loop on stator control points 
DO 10 I = 1, NP 

21 = 0 . 5D0* ( 1D0 -COS (PI* (2 DO* I - 1 DO ) / ( 2D0*NP) ) ) 
Ell = { ALF1 *CT2 + BET*ST2 ) *C2 *ZI 
E2I = { ALF2*CT2+BET*ST2 ) *C2*ZI 
E3I = ( ALF3 *CT24BET*ST2 } *C2 *ZI 

c.. Effect of rotor on stator, eq . 51 


IF 

(K .NE. 0) THEN 





KRS 

= 1/MR2 * ( 





> 

(BET*CT2-ALF1*ST2) * 

VR1 

*EXP(AI* 

( ALF1* 

(XS-XE)+E1I) ) 

> 4 

( BET*CT2 -ALF2* ST2 ) * 

(KR24VR2) 

* EX F* ( AI* 

{ ALF2* 

XS 4E2I)) 

> 4 

(ALF3*CT2+BET*GT 2) * 

(KR34VR3) 

* EXP ( AI* 

( ALF3* 

XS +E3I)) 


c.. Place elements in KMATRX , forming real elements, from complex KRS, eq. 69 
KMATRX (NPNH24 ( 2 *ABS ( K ) - 2 ) *NP+I , ~( 2 *N-2 ) *NP+J ) = 

>KMATRX(NPNH2+ (2*ABS (K) -2) *NP+I , (2*N-2 ) *NP+J ) +REAL (KRS) 

KMATRX (NPNH2+ (2*ABS(K) -2) *NP+I, (2*N-1 ) *NP+J) = 

>KMATRX (NPNH2+ { 2 *ABS { K ) -2 ) *NP+I , (2*N-1 } *NP+ J ) - I MAG ( KRS ) 

KMATRX (NPNH2+ (2*ABS{K) -1) *NP+I, ( 2*N-2 ) *NP+J ) = 

>KMATRX (NPNH2+ ( 2 *ABS (K ) - 1 ) *NP+I , ( 2 * N- 2 ) *NP+J } + 1 SIGN ( 1 , K) *IMAG ( KRS) 

KMATRX (NPNH2+ (2*ABS<K)-1) *NP+I, (2*N-1 ) *NP+J) = 

>KMATRX (NPNH2+ (2*ABS (K ) -1 ) *NP+I, { 2 *N- 1 ) *NP+J ) + 1 SIGN ( 1 , K) *REAL ( KRS) 

ENDIF 


c.. Compute the portion of the KRR coefs caused by the reflected waves, eq . 46. 
KRR = 1/MR1* 

> ( ( BET*CT1 -ALF1 *ST1 ) *VR1*EXP{AI* <ALF1*CT1 + BET*ST1 ) *C1*ZI) *EXPE1 

> + <BET*CT1-ALF2*ST1 ) *VR2*EXP{AI* { ALF2 *CT1 4BET*ST1 } *C1*ZI } 

> + ( ALF3 *CTl4BET*STl ) *VR3*EXP(AI* (ALF3*CT1+BET*ST1 } *C1*ZI ) ) 


c.. Form real elements, do sum over K, and place in KMATRX. 

KMATRX ( (2*N-2) *NP4l, ( 2*N-2 ) *NP+J ) = 

> KMATRX((2*N-2)*NP+I, (2*N-2)*NP+J) 4 REAL (KRR) ! eq. 62 

KMATRX ( ( 2*N- 1 ) *NP4l, (2*N-2 ) *NP+J ) = 

> KMATRX ( <2*N-1) *NP4l, (2*N-2) *Np4j) 4 IMAG(KRR) ! eq . 63 

10 CONTINUE 


c . . 


c . . 


20 


Fill in the remaining sections of the rotor-on- rotor quarter of the matrix 
from the second parts of Eqs. 62 and 63. 

DO 20 N = 1, NH 
DO 20 J = 1, NP 
DO 20 I = 1, NP 

KMATRX ( { 2*N- 1 ) *NP4l, (2*N-1) *NP+J) = 

> KMATRX { (2*N-2) *NP+I, (2*N-2 ) *Np4d) 

KMATRX ( ( 2*N-2 ) *NP4l, ( 2 *N - 1 ) *Np4j ) = 

> -KMATRX ( ( 2 *N- 1 ) *Np4l, (2*N-2) *NP+J) 

CONTINUE 


RETURN 

END 

c 

C 


Figure 4. (continued) 
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SUBROUTINE GETVS ( BETA , NLAM, MX , ST, CT, VI , V2 , V3 ) 
c.. Generates Vi and V2 (Smith's vl ' /beta and v2'/beta) and V3 (Smith's v3 ' / 
c.. alphaH) for either rotor waves or stator waves. 

c.. For stator waves, call with BETA (N, K) , NLAM=N*LAM2=N*B1* (MY1 -MY2 ) /MR2 , MXB, 
c.. ST=SIN (THETA2 ) , CT=COS (THETA2 ) . 

c.. For rotor waves, call with BETA { K, N} , NLAM=N*LAM1=N*B2* (MY1-MY2 ) /MR1 , MXB, 
c.. ST=SIN (THETA1 ) , CT=COS (THETA1 ) . Derivation given in appendix C. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION NLAM, MX 
DOUBLE COMPLEX ROOT, G, VI, V2 
INTEGER BETA 

ABAR = NLAM* *2 ♦ BETA* *2 + 2 . 0 DO* NLAM* BETA* ST ! eq . C-18 

E = BETA* * 2 - ABAR* (MX/CT) **2 ! eq . C-19 

c.. E < 0 for propagation, E > 0 for decay. Any E = 0 (resonance) cases 
c.. will be caught by the prior call to subroutine ALFBET. 

IF (E .LT. 0.0D0) THEN 

ROOT = CMPLX (SQRT(-E) , 0.0D0) 

ELSE 

ROOT = CMPLX (0.0D0, -SQRT ( E) ) 

END IF 

F = BETA + NLAM* ST 
G = NLAM * BETA * CT / ROOT 

VI = ( -F + G) / (2 . ODD* ABAR) ! eq . C- 15 

V2 = ( F + G) / (2 . ODD* A BAR) ! eq . C-16 

V3 = - NLAM*CT/ ABAR ! eq . C-17 

RETURN 

END 

C 

c 


Figure 4. (continued) 
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SUBROUTINE ECOEFS (NH, NR , Cl , C2 , SC2 , CT1 # CT2 , ST1 , ST2 , XS, 

> MXB, LAM2,ALF, BETA, R12 , R21 , R13 , R31 , T14 , T28 , T3 8 , KMATRX , KSUP, KSDN) 

c.. Generates elements of the matrix of influence coefficients KMATRX that give 
c . , the upwash caused by stator loading at control points on the rotor and 

c.. stator. These are computed from KSR ( N, K , I , d ) , effect of stator on rotor, 

c.. and KSS ' ( N, I , d ) 1 , effect of stator on stator. The prime on KSS' indicates 

c.. that only the part of KSS associated with waves reflected from the actuator 
c.. disk is computed here. The remainder is computed in GENKSS using original 
c.. routines from the Smith code. N counts the stator loading harmonics, I the 
c.. control points, and J the load elements. K counts the cascade wave index in 
c.. the stator frame which becomes the time harmonic index in the rotor frame. 

IMPLICIT DOUBLE PRECISION (A-H,OZ) 

DOUBLE PRECISION MXB, MR1 , MR2 , LAM2 , KMATRX ( 1 020 , 1 02 0 ) 

INTEGER BET, BETA (-5:5, - f» : r . ) 

DOUBLE COMPLEX A I , ALF1 , ALF2 , ALF3 , EXPE1 , El I , E2 1 , El I , 

> KS1 , KS2 , KS3 , VI , V2 , VS1 , VS2 , VS 3 , KSR, KSS , 

> ALF(9,-5:5,-T. : S ) , 

> KSUP ( 5 , - S : S , r . 1 ) , KSDN ( 5 , - 5 : r » , 51 ) , 

> R12 (-5 : 5, -S : 5) , R21 {-S: 5, -5 :5) , R13 (-5 : 5, -5: 5) ,R31 (-5:5, -5:5) , 

> T14 (-5 : c ), -5 : S) ,T2R (-S : S, -5 :5) , T38 (-5 : 5 , -5 : 5) 

AI = (0.0D0, 1.0D0) 

PI = 3. 141S926SD0 

WRITE ( * , *) 'Entering 
MR1 = MXB/CT1 

MR2 = MXB/CT2 

CON = MR2/SC2 

XE = XS + C2*CT2 

NPNH2 = NP* NH * 2 
DO 1 0 N = 1 , NH 

DO 10 K = -NH, NH 
BET = BETA ( N , K ) 

ALF1 = ALF ( 1 , N, K ) 

ALF2 = ALF { 2 , N, K } 

ALF 3 = ALF ( 3 , N, K) 

CALL GETVS ( BET , N* LAM2 , MXB , ST2 , CT2 , VI , V2 , V3 ) 

DO 10 J = 1, NP ! Loop on load elements 

ZOd = 0.5DO*(1DO - COS ( PI* ( d-lDO) /NP) ) ! Chordwise locations for loads 

KS1 = CON* VI *EXP ( - AI * ( ALF1 *CT2 + BET*ET2 ) *C2 *Z0d ) ! 

KS2 = CON*V2*EXP(-AI* (ALF2*CT2+BET*ST2) *C2*Z0d) ! eq. 22 

KS3 = CON*V3*EXP(-AI* { ALF3 *CT2+BET*ST2 ) *C2*Z0d) ! 

c.. Compute VS1 , VS2, VS3 from eq , 26 

c.. Hold out EXPEl=exp { -AI *ALFI *XE) from VS1 here to avoid overflow later in KSR 
EXPE1 = EXP ( -AI *ALF1*XE} 

VSI = ( ( R12 (N, K } *R21 (N,K) *EXP ( AI*ALF2*XE) 

> 4 R13(N,K)*R31 (N, K) * EXP ( AI *ALF3 *XE ) ) *KS1 *EXP { - AI *ALF1 *XS ) 

> 4 R21 (N, K) *KS2*EXP(AI*ALF2* (XE-XS) ) 

> 4 RU (N,K) *KS1*EXP(AI*ALF3* (XE-XS) ) ) 

> / (1.0 DO - R12(N,K)*R21(N,K)*EXP<AI*(ALF2-ALF1}*XE) 

> - R1 UN, K) *R 31 (N, K ) *EXP ( AI * ( ALF3-ALF1 ) *XE) ) 

VS2 = ( KS1 * EXP ( - AI *ALF! *XS) +VS1*EXPE1 ) * R1 2 ( N, K ) 

VS 3 = ( KS1 * EX P { - AI * ALF1 *XS) +VS1*EXPE1 ) *R13 (N, K) 

c.. KSUP and KSDN are passed out <»f the subroutine for later use in computing 
c.. pressure in the non-swirl regions upstream and downstream. 

KSUP {N, K, d ) = (KS)*EXP(-AI*ALF1*XS) +VS1*EXPE1) * T14(N,K) ! eq . 77 

KSDN (N, K, d ) = (KS2*EXP(-AI*ALF2*XS)+VS2)*EXP(AI*ALF2*XE)*T28(N,K) ! eq . 88 
> 4 ( KS 3 * EXP ( -AI * ALF3 *XS } +VS3 ) * EX P { AI *ALF3*XE) *T3 8 (N, K) 


! Relative Mach #, row 1 
! Relative Mach #, row 2 
! A constant 

! Axial coordinate of stator exit 
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c.. Loop on rotor control points 
DO 10 I = 1, NP 

ZI = 0.5D0* (1DO-COS(PI* (2*I-1D0) / (2*NP) ) ) 

Ell = ( ALF1 *CT1 +BET*ST1 ) *C1 *ZI 
E2I = (ALF2*CT1+BET*ST1 ) *C1 *ZI 
E3I = (ALF3*CTl4BET*STl)*Cl*ZI 

c . . Effect of stator on rotor, eq . 37 
IF (K .NE. 0) THEN 
KSK = 1/MR1* ( 

> (BET *CT1 - ALF1 *ST1 ) *KS1 *EXP ( AI * El l ) *EXP ( -AI *ALF1 *XS ) 

> + (BET *CT1 - ALF1*ST1) *VS1 *EXP (AI*E1 I ) *EXPE1 

> + (BET *CT1 - ALF2*ST1) *VS2*EXP(AI*E2I) 

> + ( ALF3 *CT1 + BET* ST1 ) *VS3*EXP (AI*E3I ) ) 

c.. Place elements in KMATRX , forming real elements, from complex KSR, eq . 70 
KMATRX ( (2*ABG(K) -2) *NP+I,NPNH2+ ( 2 *N- 2 ) *NP4J) = 

> KMATRX ( (2*ABS (K) -2 ) *NP+I , NPNH2+ ( 2*N-2 ) *NP4J) +REAL(KSR) 

KMATRX ( ( 2 *ABS (K ) -2) *NP+I,NPNH2+ (2*N-1 ) *Np4j) = 

> KMATRX ( ( 2 * ABK ( K ) -2 ) *Np4l , NPNH2+ ( 2*N- 1 ) *NP4 J ) -IMAG (KSR) 

KMATRX ( ( 2 *AB! : ! (K ) - 1 ) *NP+I . NPNH2+ (2*N-2 ) *Np4j) = 

>KMATRX ( ( 2 *AB.' : ! (K ) -1) *Np4 1 , NPNH2 + ( 2 *N-2 ) *NP+J ) +ISIGN ( 1 , K) * IMAG (KSR) 

KMATRX ( (2*AB:; (K) -1 ) *Np4l, NPNH24 (2*N-1 ) *NP+J) = 

> KMATRX ( (2*AB: : ; (K) -1 ) *Nt‘4l , NPNH24 ( 2*N- 1 ) *NP+J) +ISIGN ( 1 , K ) *REAL ( KER ) 

ENDIF 

c.. Compute the portion of the KSS coefs caused by the reflected waves, eq. 31. 
KSS = 1/MR2* 

>( (BET*CT2 - ALF1 * ST2 ) *VS1 * 

> EXP ( AI * ( ALF1 * (XS-XE)4(ALF1*CT24BET*ST2)*C2*ZI) ) 

> 4 ( BET*CT2 - ALF2 * ST 2 ) *VS2 * 

> EXP (AI * (ALF2*XS4 (ALF2*CT24BET*ST2) *C2*ZI) ) 

> 4 (ALF3*CT2 4 BET*GT2) *VS3* 

> EXP ( AI * 1 ALF < *XS4 (ALF3*CT24BET*ST2) *C2*ZI ) ) ) 

c. . Form real elements, do sum over K, and place in KMATRX 
KMATRX (NPNH24 (2*N-2 ) *Np4l , NPNH24 ( 2*N-2 ) *NP4J) = 

> KM/.TRX (NPNH24 (2*N-2) *NE'4l, NPNH24 (2*N-2) *NP4d) 4 REAL (KSS ) ! eq . 64 

KMATRX (NPNH24 (2*N-1) *NP4l ,NPNH24 (2*N-2) *NP+J) = 

> KMATRX(NPNH24(2*N-1)*NP4I,NPNH24(2*N-2)*NP4J) + IMAG(KSS) ! eq. 65 

10 CONTINUE 

c.. Fill in the remaining sections of the stator on stator quarter of the matrix 
c.. from the second parts of Eqs. 64 and 65. 

DO 20 N = 1, NH 
DO 20 J = 1, NP 
DO 20 I = 1, NP 

KMATRX (NPNH2+ (2*N-1 ) *Np4l , NPNH2+ ( 2*N-1 ) *Np4d ) = 

> KMATRX (NPNH2 4 ( 2*N-2 ) *Np4 I , NPNH24 (2*N -2 ) *NP+J ) 

KMATRX (NPNH24 (2*N-2) *Np4l, NPNH24 (2*N-1 ) *Np4.1) = 

> -KMATRX (NF'NH2 4 ( 2*N- 1 ) *NF'4 I , NPNH24 (2*N-2 ) *Nf'4j ) 

20 CONTINUE 

RETURN 

END 


Figure 4. (continued) Source code for CUP2D. 
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SUBROUTINE GENKRR (NH , NP , Cl , C2 , SCI , SC2 , CT1 , ST1 , MXB , LAM1 , KMATRX ) 
c.. This generates the input needed to call Smith's matrix generation routines 
c.. for the effect of the rotor on itself via the direct waves. It fills the 
c.. WHEAD common block and calls Smith's routine DSWK . This returns the real 
c.. and imaginary parts of the matrix. These are placed in the appropriate 
c, , locations in KMATRX, adding to the elements already computed by RCOEFS that 
c.. account for the effect of the rotor on itself via the reflected waves. 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /WHEAD/ SC, STAG, MACH, LAM, PHASE, DEG, PI , COSST, SINST, 

> MACH2 , B, BETA 2 , BC ( BC2 

DOUBLE PRECISION MXB, LAM1, KR { 51 , 51 ) , KI ( 51 , 51 ) , KMATRX ( 1 02 0 , 1 02 0 ) 
DOUBLE PRECISION MACH , LAM , MACH2 


c . . 


WRITE { 

★ 

, *) 'Entering GENKRR' 

Fill /WHEAD/ common block (except for 

SC 

= 

SCI 

STAG 

= 

ASIN(STl) 

MACH 

= 

MXB/CT1 

MACH2 

= 

MACH* * 2 

BETA2 

=: 

1D0 - MACH 2 

B 

- 

SORT ( BETA2 ) 

DEG 

= 

57 .29578D0 

PI 


3 .14159265D0 

COS ST 

= 

CT1 

SINST 

= 

ST1 

BC2 

= 

1D0 - MXB* *2 

BC 

= 

SQRT ( BC2 ) 

Loop) on 

rotor loading harmonic: 

DO 300 

N=l, NH 

LAM = 

N 

*LAM1*C1 

PHASE 

= 

2 . 0 D 0 * F ’ I * N * C 1 / C 2 * S C 1 / S C 2 


PHASE) 


c.. Call Smith's original matrix generation routine 
CALL DSWK (KR, KI , NP, IW) 

IF (IW . EQ . 1) THEN ! Smith's resonance check 

WRITE ( * , 1 ) N 

1 FORMAT {IX,' DSWK RETURNED IW=1 TO GENKRR FOR N=',I2) 

STOP 'EXECUTION TERMINATED DUE TO RESONANCE' 

ENDIF 


c . . Add 


> 


> 


100 

300 


Smith's real and imaginary matrix elements 
DO 100 I = 1, NP 
DO 100 J = i, NP 

KMATRX ( ( 2 *N-2 ) *NP+I , ( 2*N-2 ) *NP+ J) = 
KMATRX ( (2*N-2) *NP+I, (2*N-2) *NP+J) + 
KMATRX ( (2*N-1 ) *NP+I , ( 2*N- 1 ) *NP+J ) = 
KMATRX ( (2*N-1 ) *NP+I, (2*N-1 ) *NP+J) + 
KMATRX ( (2*N-1 } *NF>+I , (2*N-2) *NP+J) = 
KMATRX ( (2*N-1 ) *NF’+ I , ( 2*N-2 ) *NP+J ) + 
KMATRX { { 2 *N-2 ) *NP+I , ( 2*N- i ) *NP+J ) = 
KMATRX ( ( 2*N- 2 ) *NP+I, (2*N-1) *NP+J) - 
CONTINUE 
CONTINUE 


into KMATRX, 


KR ( I , J ) 
KR ( I , J ) 
KI (I, J) 
KI (I, J) 


see Eqs . 62 & 63 


RETURN 

END 


Figure 4. (continued) Source code for CUP2D. 
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SUBROUTINE GENK SS ( NH , NI> , C 1 , C2 , SCI , SC2 , CT2 , ST2 , MXB , LAM 2 , KMATRX ) 

c.. This generates the input needed to call Smith's matrix generation routines 
c.. for the effect of the stator on itself via the direct waves. It fills the 
c.. WHEAT common block and calls Smith's routine DSWK . This returns the real 
c.. and imaginary parts of the matrix. These are placed in the appropriate 
c.. locations in KMATRX, adding to the elements already computed by SCOEFFS that 
c.. account for the effect of the stator on itself via the reflected waves. 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /WHEAT/ SC, STAG, MACH , LAM, PHASE, TEG, PI , COSST, SINST, 

> MACH2 , B, BETA2, BC, BC2 

DOUBLE PRECISION MXB, LAM 2 , KR ( 51 , 51 ) , KI ( 51 , 51 ) , KMATRX (1020, 1020) 

DOUBLE PRECISION MACH , LAM , MACH2 

WRITE (*,*)' Entering GENKSS ' 

C. . . Fill /WHEAT/ common block (except for LAM and PHASE) 


SC 

= 

SC2 

STAG 

— 

ASIN ( ST2 ) 

MACH 

= 

MXB/CT2 

MACH2 

= 

MACH* * 2 

BETA2 

= 

1.0 DO - MACH 2 

B 


SQRT ( BETA2 ) 

DEG 

= 

57 . 2*V78D0 

PI 

- 

3 .1415^26500 

COSST 


CT2 

SINST 

- 

ST2 

BC2 

= 

1 .0D0 - MXB** 

BC 

= 

SQRT ( BC2 ) 


C... Loop on stator loading harmonic 
NPNH2 = NP*NH*2 
DO 300 N= 1 , NH 
LAM = N* LAM2 *C2 
PHASE = -2 . 0D0* PI*N*C2/C1 *SC2/SC1 

C... Call Smith's original matrix generation routine 
CALL DSWK (KR, KI~, NP, IW) 

IF (IW . EQ. 1) THEN ! Smith's resonance check 

WRITE ( * , I ) N 

1 FORMAT ( IX , ' DSWK RETURNED IW=1 TO GENKSS FOR N=',I2) 

STOP 'EXECUTION TERMINATED DUE TO RESONANCE' 

ENTIF 

C. . Add Smith's real and imaginary matrix elements into KMATRX, see Eqs . 64 & 65 
TO 100 I = 1, NP 
TO 100 J = 1, NP 

KMATRX (NPNH2+ ( 2 *N- 2 ) *NP+I ,NPNH2+ (2*N-2) *NP+J) = 

> KMATRX (NPNH2+ ( 2*N-2 ) *NP+I , NPNH2+ (2 *N-2 ) *NP+J ) + KR(I,J) 

KMATRX (NPNH2+ (2*N-1 ) *NP+I , NPNH2+ ( 2 * N- 1 ) *NP+J) = 

> KMATRX (NPNH2+ I 2*N-1 1 *NP+I , NPNH2+ (2*N-1 ) *NP+J) + KR(I,J) 

KMATRX (NPNH2+ ( 2*N- 1 1 *NP+I ,NPNH2t (2*N-2) *NP+J) = 

> KMATRX (NPNH2+ ( 2*N-1 ) *NP+I , NPNH2+ ( 2 * N - 2 ) *NP+J) + KI(I,d) 

KMATRX (NPNH2+ (2*N-2) *NP+I ,NPNH2+ ( 2 * N- 1 ) *NP+,I) = 

> KMATRX (NPNH2+ ( 2 * N - 2 ) *NP+I , NPNH2+ { 2 * N - 1 ) *NP+J ) - KI(I,J) 

100 CONTINUE 

300 CONTINUE 

RETURN 

ENT 

C 

c 


Figure 4. (continued) Source code for CUP2D. 
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SUBROUTINE SOLVE(NH,NP, KMATRX , WR , WS f LR , LS ) 

c.. This-, routine finds the loading on both the rotor and the stator 

c.. simultaneously by matrix inversion. Externally imposed upwash velocities 

c.. enter the routine in complex form (WR,WS) and are used to form a one 

c.. dimensional real vector WASH. KMATRX is inverted using the LINPACK routines 

c.. and multiplied by WASH. The result is the one dimensional load vector LOAD. 

c.. This is decomposed into the complex load vectors LR and LS and these are 

c,. sent to the LOADS routine for output. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION KMATRX (1020,1020) , LOAD(1020), WASH(1020) 

DOUBLE COMPLEX WR(5,51>, WS{5,51), LR(5,51), LS (5,51) 

ZERO = 0 . 0D0 

NPNH2 = NP*NH*2 ! constant needed for element shifting 

NPNH4 = NP*NH*4 ! size of one side of real matrix 

c.. Invert matrix of influence coefficients 
CALL MATINV (KMATRX, NPNH4 ) 
c.. KMATRX is now KMATRX inverse 

c.. Form real upwash vector from complex upwash vectors for rotor and stator, 
c.. Minus signs because objective is to find the loading that produces upwash 
c.. to cancel wake upwash. Note, rotor upwash is zero for the usual rotor/ 
c.. stator interaction problem, but the code could deal with rotor vibration or 
c.. stator potential field excitation also. 

DO 10 K = 1, NH 
DO 10 I = 1, NP 

WASH ( { 2*K-2 ) *NP+ I ) = - REAL (WR(K,I)} ! eq . 71 

WASH ( { 2 * K - 1 ) *NP+I) = - IMAG (WR(K,I)) ! 

WASH ( NPNH2 + ( 2*K- 2 ) *NP+ I ) = - REAL (WS{K,I)) ! 

WASH { NPNH2 + ( 2 * K - 1 ) *NP+ 1 ) = - IMAG (WS(K,I)) ! 

10 CONTINUE 

c.. Zero load vector before matrix multiplication 
DO 20 NU = 1, NPNH4 
LOAD (NU) = ZERO 
20 CONTINUE 

c.. Matrix multiplication to get load vector, which is in real form. 

DO 30 NU = 1, NPNH4 
DO 30 MU = 1, NPNH4 

LOAD { NU ) = LOAD ( NU ) + KMATRX ( NU , MU ) * WASH (MU) 

30 CONTINUE 

c.. Form complex load vectors for rotor and stator from single real load vector 
DO 4 0 N = 1 , NH 
DO 4 0 .1 =: 1, NP 

LR ( N , i7 ) =CMPLX ( LOAD ( ( 2*N- 2 ) *NP + J ) , LOAD ( ( 2*N- 1 ) *NP+J ) ) ! eq . 

LS (N, J) =CMPLX (LOAD(Nf , NH24 (2*N- 2 ) *NP + J ) , LOAD (NPNH2+ (2*N-1) *NP+J) ) ! 72 
40 CONTINUE 

CALL LOADS ( NH , N P , LR , LS ) 

RETURN 

END 


Figure 4. (continued) Source code for CUP2D. 
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SUBROUTINE LOADS { NH , NP, LR, LS) 

c.. Computes and prints integrated loading on both blade rows for multiple 
c.. harmonics ♦ Formulas equivalent to Smith's eq . 56. Values are lift per 
c.. unit span divided by chord* rho* relv* *2 . Treatment of first and last terms 
c.. handled automatically by inversion process. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE COMPLEX LR(5,51), LS (5,51) , CLR, CLS 

WRITE ( * , * ) ' ' 

WRITE ( * , *} ' 

WRITE { * , * ) ' N 

WRITE ( * , * ) ' 

DO 20 N=l, NH 

CLR = (0D0, 0D0) 

CLS = (0D0, 0D0) 

DO 10 J = l, NP 

CLR = CLR - LR ( N , J ) 

CLS = CLS - LS(N,J) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
c.. Turn on the next 2 lines for printout of chordwise load distributions 
C WRITE ( * , 5 ) N, 0, LR ( N , U ) , LS(N,J) 

C5 FORMAT ( IX , 12,15,2015.5, ' ' ,2015.5) 

ccccccccccccc:ccc:ccccccccccccccc:ccccccc:cr:c:cccccccccccccccccccccccccccccccc 

10 CONTINUE 

WRITE {*,!) N, CLR, CLS 

1 FORMAT (IX, 12, ' \ 2F10.5,' \2F10.5 ) 

20 CONTINUE 

WRITE (*,*} ' ' 

RETURN 
END 

C 

C 


f * VALUES OF 
CLROTOR (N) 
real imaq 


LIFT COEFFICIENTS ***' 

CLSTATOR ( N ) ' 
real imag 


Figure 4. (continued) 
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SUBROUTINE OUTPUT (NH, NP , B1 , B2 , MXA , MXB , MXC, MS , MY1 , MY2 , ALF, ABA, ABC, 

> XS,XA, XC, POPSA, POPSC,KRUP, KRDN, KSUP, KSDN, LR, LS ) 

c.. Calculates sound pressure at axial locations XA and XB and sound power 
c.. (average per unit area) upstream and downstream based on loading from 
c.. SOLVE and influence functions from INFFNS. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION MXA, MXB, MXC, MS, MY1 , MY2 , PWLAT ( 5 ) , PWLCT ( 5 ) 
INTEGER Bl, B 2, ORDER 

DOUBLE COMPLEX A4,A8,G4,G8, PUP, PDN , ALF (9,- 5:5, -5:5) , AI , 

> LR (5, SI), LS (5,51), LUP (5, -5:5), LDN (5, -5:5), 

> LRUP ( 5 , - 5 : 5 ) , LRDN ( 5 , - 5 : 5 ) , LSUP { 5 , - 5 : 5 ) , LSDN ( 5 , - 5 : 5 ) , 

> KRUP (5, -5 : 5 , SI) , KRDN (5, -5: 5, 51 ) , KSUP (5, -5:5,51), KSDN (5, - 5 : 5, 51 ) 

AI = ( 0 . 0D0 , 1 . 0D0) 

WRITE (*,*) ' ' 

WRITE (*,*) 'Axial locations for sound pressure output in radii from 

> rotor leading edge ' 

WRITE ( * , 5 ) XA~ XC 

5 FORMAT (IX, ' F» u PRESup, Xa = ' , F5 . 3 , ' ( For PRESdn, Xc = ' , F5 . 3 ) 

WRITE ( * , *) ' ' 

WRITE (*, * )' Decibel Levels for Pressure Waves and Power Levels 

> I Cutoff Ratios' 

WRITE (*,*)' N K FREQ nBl-kB2 PRESup PRESdn PWLup PWL 
>dn I A B C' 

c.. Zero the wave accumulators 
DO 10 N = 1, NH 

DO 10 K r -NH, NH 

LRUP (N, K) = (0D0, 0D0) 

LRDN ( N , K ) = (0D0, 0D0) 

LSUP (N, K ) = (ODD, 0D0) 

LSDN ( N, K ) = (0D0, 0D0) 

10 CONTINUE 

c.. Sum (Contributions to waves over load elements 
DO 20 N = 1, NH 

DO 20 K = -NH, NH 

DO 20 J = 1, NP 

LRUP (N , K ) = LRUP (N, K ) + KRUP {N, K , J ) *LR (N , J ) ! eq . 84 

LRDN (N, K ) = LRDN (N, K ) + KRDN (N, K , d ) *LR (N, J ) ! eq. 96 

LSUP { N , K ) = LSUP (N , K ) + KSUP ( N , K , J ) * LS ( N , J ) ! eq . 77 

LSDN (N , K ) = LSDN ( N , K ) + KSDN (N, K , J ) *LS (N , J ) ! eq . 89 

20 CONTINUE 

c.. Add rotor waves to stator waves upstream and downstream 
DO 30 N = i , NH 
DO 32 K = 1 , NH 

LUP ( N , K) = LRUP ( K , N) + LSUP(N, K) ! eq . 106 

LUP ( N , -K ) r -CONdG (LRUP (K, - N) ) + LNUP(N,-K) 

LDN ( N , K) = LRDN ( K , N) + LSDN ( N , K) ! eq. 116 

LDN ( N , - K ) = -CONdG (LRDN (K, -N) ) + LSDN(N,-K) 

32 CONTINUE 

c.. Following terms computed without steady loadinq effect from rotor 
LU P ( N , 0 ) = LSUP ( N, 0 } 

LDN { N , 0) = LSDN (N, 0) 

30 CONTINUE 

c.. Compute modal and total powers upstream and downstream 
DO 50 N = 1, NH 

PWRAT = 0.0 [)0 ! Total power accumulator, upstream 

PWROT = O.ODO ! Total power accumulator, downstream 


Figure 4. (continued) Source code for CUP2D. 
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DO 52 K = -NH, NH 

FREQ = N* B1 *MY1 - K*B2*MY2 

ORDER = N* Bi - K*B2 ! Negative of beta(n,k) 


c.. CTRAT's below are cutoff ratios (infinity is rigged to print 9.99), eq . B-6 


IF (ORDER 
CTRATA 
CTRATB 
CTRATC 
ELSE 

CTRATA = 
CTRATB = 
CTRATC = 
CTRATA = 
CTRATB = 
CTRATC = 
ENDIF 


EQ. 0) THEN 
= 9.99D0 
= 9.99D0 
= 9.99D0 

SQRT ( (ABA* FREQ) **2/ ( ( 1D0-MXA**2) *ORDER**2) ) 

SQRT ( ( FREQ-ORDER* MS ) **2/ ( (1D0-MXB**2) *ORDER**2) 
SQRT ( ( ABC*FREQ) **2/ ( ( 1D0-MXC* *2 ) *ORDER* *2 ) ) 

MIN (ABS (CTRATA) ,9.99D0) 

MIN (ABS (CTRATB) , 9.99D0) 

MIN (ABS (CTRATC) , 9 . 99D0) 


c.. Output only waves with cutoff ratios > 0.2 in up or downstream region 
IF (CTRATA .GT. 0.2D0) GOTO 51 
IF (CTRATC .LE. 0.2D0) GOTO 52 


51 A4 = ALF (4 , N, K) 

A8 = ALF ( 8 , N, K ) 

G4 = - (ABA*FREQ + MXA*A4 ) 
G8 = - (ABC* FREQ + MXC*AH) 


! axial wavenumber upstream 
! axial wavenumber downstream 
! eq. 81 
! eq. 93 


CALL GETPWL (CTRATA, MX A, A4 , G4 , LUP ( N, K ) , - 1 , PWRA, PWLA) 
CALL GETPWL ( CTRATC, MXC, A8,G8, LDN(N, K) , 1, PWRC, PWLC) 
PWRAT = PWRAT + PWRA 
PWRCT = PWRCT * PWRC 


c.. Compute pressures at phi=0 and x = xa and xc 
PUP = G4*LUP (N, K ) *EXP { AI *A4*XA) 

PDN = G8*LDN (N, K) *EXP(AI*A8* (XC-XE) ) 

WRITE (*,i) N, K, FREQ, ORDER, DBEL ( PUP, POPSA) , DBEL ( PDN, POPSC) , 
> PWLA, PWLC, CTRATA, CTRATB, CTRATC 

1 FORMAT ( IX , 215, F7.2, 16, 4F8.1, 1 ', 3F6.2) 

52 CONTINUE 

PWLAT(N) = MAX { 0D0 , 10D0*LCG10 ( PWRAT*1 . 0D13+1 . 0D-30 ) ) 

PWLCT(N) = MAXtODO, 10D0* LOGIC) ( PWRCT* 1 . 0D13+1 . 0D-30 ) ) 

WRITE ( * , 3 ) N , PWLAT ( N) , PWLCT(N) 

3 FORMAT ( IX , ' Total power for N=',I2,' \2F8.1) 

WRITE { * , * ) ' ' 

50 CONTINUE 

RETURN 

END 


Figure 4. (continued) Source code for CUP2D. 
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SUBROUTINE CETPWL (CTRAT, MX, A, G, L, IX, ROWER, PWL) 
c.. Computes sound power level according to theory in Eqs . 98-116. 
c.. Checks that power is real (within numerical accuracy) and that real part 
c.. has the correct sign for flux away from blades. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION MX 
DOUBLE COMPLEX A, G, L, PWR 

c.. No power for waves that are out off. 

IF (CTRAT . LE. IDO) THEN 
POWER = ODD 
PWL = 0[)0 
RETURN 
ENDIF 

c.. Treat power as complex to verify correct behavior. See Eqs . 112 & 115. 
c.. SIGN function below is needed because derivation was for power flux in the 
c.. +x direction. The expression below must be real and > 0 for regions A & C. 

PWR = ( (1D0+MX**2) *A*G a MX* ( A* *2 +G* * 2 ) ) *ABS { L ) **2 *SIGN(1,IX) 
c.. Check that PWR is real. 

TEST = ABS ( IMAG (PWR)/ (REAL (PWR) +1 .0-20} ) 

IF (TEST . GT . . OOIDO) THEN 

STOP 'Execution terminated because power flux is complex' 

ENDIF 

c.. Check that power flux is outgoing on either side of the source. 

IF (REAL ( PWR } . LT. 0.) THEN 

WRITE { * , * ) ' IX - ' , IX 

STOP 'Execution terminated because power flux is negative' 

ENDIF 

POWER = REAL (PWR) 

IF (POWER .LT. 1.0D-13) THEN 

PWL = 0D0 

ELSE 

PWL = 10 . 0D0*LO(.;i0 ( POWER* 1 . ODi A ) ! eq . 113 

ENDIF 

RETURN 

END 

C 


DOUBLE PRECISION FUNCTION DBEL(X,POPS) 

c.. Computes SPL dB from harmonic peak value normalized by ambient rho*c:**2. 
c.. Uses POPS, the ratio of local ambient pressure to 2116 psf. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE COMPLEX X 

TEMP=20 . 0 DO* LOG 10 ( . 7071 1D0* ABS ( X ) *i .400 * POPs’* 5 . 0651D9+1 .D-35) 

DBEL=MAX (TEMP , 0 . 0D0 ) 

END 


Figure 4. (continued) Source code for CUP2D. 



C Calculation of Kernel Matrix Elements from Original Smith Code 
C 

SUBROUTINE DSWK (KR, KI , NP, IW) 

C 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION MACH , LAM , MACH2 , MACH4 , MACH6 , KR ( 51 , 51 ) ,KI (51, 51} 
COMMON /WHEAD/ SC , STAC , MACH , LAM , PHASE , DEG , PI , COSST , SINST , 

> MACH2, B,B2,BC, BC2 

COMMON / WAVEC / XU # APU , API), ANU, AND, PUR, PDR, PUI, VU,VD 
DIMENSION I CHECK (51,51) ,ZE(51) , ZP{51) 

C 

C.... CONSTANTS FOR VORTEX SHEET CALCULATION 
X = LAM* SC* COS ST 

Y = LAM*SC*SINST + PHASE 

VORT = 0 . 5* LAM* SINH (X } / (COSH { X) - COS(Y)) 


C 

C. 


C 

c. 


c 

c. 


c 

c 

c 

c 


CONSTANTS FOR LOG SINGULARITY CORRECTION 
MACH4 = MACH 2 * MACH 2 
MACH 6 = MACH2 *MACH4 

B4 = B2 * B2 

B6 = B2 * B4 

A1 = 1.0 - 0 .5*MACH2/B2 

A2 = 1.0 - 0.5/B2 4 0.9*MACH2/B4 

A3 = 0.5* (1.0 - 1.0/B2 4 MACH2 / ( 6 . 0*B4 ) 4 1.0/(3.0*B4) 
- 0.375* MACH4 / B 6 4 MACH6/ ( 6 . 0*B6) ) 


. MATCHING AND VORTEX POINTS 
DO I = 1,NP 

EPSIL = PI* FLOAT (2*1 - 1 ) /FLOAT { 2*NP) 
ZE ( I ) = 0.5* (1.0 - COS (EPSIL}) 

PSI = PI* FLOAT (I - 1) /FLOAT (NP) 

ZP(I) = 0.5* (1 .0 * COS (PSI) ) 

ENDDO 


ZERO COUNTS AND ARRAYS 
IR = 0 
I COUNT = 0 
IW = 0 
NP2 = NP*NP 
DO I = 1,NP 
DO J = 1 , NP 

I CHECK (1,0) = 0 
KR(I,J) = 0.0 
K I (1,0) = 0.0 
ENDDO 

endix:) 

ASSEMBLE MATRIX 

I ( = M 4 1 IN PAPER) GIVES VORTEX POSITION 
J( = L 4 1 IN PAPER ) GIVES MATCHING POINT 


c 

30 


C 

c 


CALL WAVE ( I R , IW) 

IF(IW.EQ.l) RETURN 
DO I = l’,NP 
DO J = 1,NP 

I F ( ICHECK ( I , J ) . EQ . 1 ) GO TO 131 
POS = ZE ( I ) - ZP(J) 

I F { POS . GT .0.0) GO TO 90 

UPSTREAM POINT 

XP = EXP ( XU* POS ) 

YP = APU* POS 
OR = XP*COG ( YP) 

QI = XP*SIN(YP) 

TERMR - { PUR*QR - PUI*QI)/SC 


Figure 4. (continued) Source code for CUP2D. 
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TERMI = ! PUR *QI + PUI*C>R)/SC 
GO TO 100 


C 

c DOWNSTREAM POINT 

90 XP = EXP ( - XU*POS) 

YP = A PD* POS 


QR = XP*COS ( YP) 

QI = XP*SIN (YP) 

TERMR = ( PDR*QR - PUI*QI)/SC 
TERMI = (PDR*0I + PUI*0R)/SC 
C 


C ADD TO MATRIX 

100 KR ( I , 0 ) = KR ( I , J ) + TERMR 

KI (1,0) = KI(I, 0) + TERMI 
C 

C CHECK CONVERGENCE OF SERIES 

C.. The next. 3 lines modified from Smith's code on 8/19/91 by DBH 

X = ABS (TERMR) + ARS (TERMI) 

Y = ABS ( KR (1,0) ) + ABS ( KI (1,0) ) 

IF( (X/Y) .GT.1.0D - 7) GO TO 131 
C X = TERMR*TERMR + TERMI*TERMI 

C Y = KR (I ,0) *KR (1,0) + KI (1,0) *KI (I ,0) 

C IF ( (X/Y) .GT.1.0D - 11) GO TO 131 


I CHECK (1,0) = 1 
I COUNT = I COUNT + 1 
C 

C CORRECT FOR LOG SINGULARITY (LAST TIME THROUGH) 

SUM = 0.0 

EPSIL = FT * FLOAT (2*1 - 1 ) /FLOAT ( 2 *NP) 

PSI = FT * FLOAT (0 - 1) /FLOAT (NP) 

NPM1 = NF* - 1 
DO OR = 1 , NF’Ml 

FOR = FLOAT (OR) 

SUM = SUM + COS (F0R*EPSIL) *COS (FJR*PSI) /FOR 
ENDDO 

SUM = 2.0* SUM + LOG (4 . 0*ABS (POS) ) 

SUM = SUM * LAM / ( 2 . 0 * PI * B ) 

PLAM = LAM* POS 
PLAM2 = PLAM* F’LAM 
PLAM 3 = F’LAMi * PLAM 

KR ( I , 0 ) = KR(I,J) + SUM* (Al* PLAM - A3*PLAM3) 

K I ( 1 , 0 ) = K I ( 1 , 0 ) + SUM* (1.0 - A2 * F’LAM2 ) 

C 

C ADD VORTICITY WAVE 

IF ( POS . LE . 0 . 0 ) GO TO 131 

KR (1,0) = KR ( I , J ) + VORT*COS ( PLAM) 

K I (1,0) = KI (I, J) - VORT*SIN ( F’LAM) 

131 CONTINUE 

ENDDO 
ENDDO 
C 

C.... CHECK FOR COMPLETION 

I F ( I COUNT . Eg . NJ’2 ) RETI JRN 
IF(IR.GT.O) THEN 
IR = -IR 
ELSE 

IR = -IR 4 1 
ENDIF 
GO TO 30 
END 


Figure 4. (continued) Source code for CUP2D. 
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c 

C CALCULATION OF ACOUSTIC WAVE PROPERTIES 
C 

SUBROUTINE WAVE(IR.IW) 

C 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION MACH , LAM , MACH 2 

COMMON /WHEAD/ SC, STAG, MACH , LAM, PHASE, DEG, PI , COSST, SINST, 
> MACH2 , B, B2 , BC, BC2 

COMMON /WAVEC/ XU, APU, APD, ANU, AND, PUR, PDR, PUI,VU,VD 
C 

BETAH=( PHASE-2.0* PI * FLOAT (IR) ) /SC 
BET AH 2 = BETAH* BETAH 

A=LAM*LAM+BETAH2+2 . 0*LAM* BETAH* SINST 

D=MACH2 * ( LAM+BETAH*SINST) *COSST/BC2 

E=BETAH2-MACH2*A 

I F ( E . NE . 0 . 0 ) GO TO 3 2 

WRITE (*,31) IR 

31 FORMAT!' Resonance at IR=',I4) 

IW=1 

GO TO 60 

32 F=SQRT ( ABS ( E) ) 

FB=F/BC2 

H= ( BETAH+ LAM * S I NST ) / ( 2 . 0 * A ) 

P=BETAH * LAM *COSST / ( F*2 . 0*A) 

IF(E.GT.O.O) GO TO 50 
C 

C WAVE NUMBERS, PROPAGATING CASE 
C 

ACUI=D+FB 

ACDI=D-FB 

A PU = ACU I * COS ST+B ETA H * S I NST 
APD=ACDI*COSST+ BETAH* S INST 
ANU = B ETAH * COS ST - ACU I * S I N ST 
AND=B ETAH * COS ST - ACD I * S I NST 
PUR = ANU* (P-H) 

PDR=AND* (P+H) 

PUI=0.0 

XU=0.0 

VU= (P-H) * (LAM+APU) /SC 
VD= (P+H) * (LAM+APD) /SC 
GO TO 60 
C 

C WAVE NUMBERS, DECAYING CASE 
C 

50 APU=D*COSST+ BETAH *SINST 

APDsAPU 

ANU = BETAH * COSST - D* S I NST 
AND=ANU 

PUR = - ANU * H - FB * S INST * P 
PDR =- PUR 

PU I = ANU * P - F B * S I NST * H 
XU= FB*COSST 
60 RETURN 

END 
C 


Figure 4. (continued) Source code for CUP2D. 
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SUBROUTINE MATINV ( A , N ) 

C. . This routine written by D.B. Hanson to call the LINPACK routines for 
C. . inversion of real matrices. Call for matrix A. Inverse returned in 
C . . same array . 

DOUBLE PRECISION A (1020, 1020) , WORK (1020) ,DET{2) , RCOND, Z ( 1 02 0 ) 
INTEGER I PVT (1020) 

WRITE (*,*) ' Entering MATINV' 

CALL DGECOfA, 1020, N, IPVT, RCOND, Z) 

WRITE {*,*) 'Condition number of KMATRX = ', RCOND 
CALL DGEDI (A, 1 020 , N , I PVT, DET, WORK , 1 ) 

RETURN 

END 


Figure 4. (continued) Source code for CUP2D. 
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Section 5 
Array Dimensions 

This section shows how the arrays are dimensioned in case they need to be changed to 
accommodate more harmonics or panels on the blades. Interpret N h and N p below to be the 
maximum permitted values of number of harmonic and number of panels. 


ALF( 9. -N h :N h , -N h :N h ) 
BETA( -N h :N h , -N h :N h ) 


KRUP, KRDN, KSUP, KSDN( N h , -N h :N h , N p ) 

R12, R21, R13, R3I, T14, T28, T38( -N h :N h , -N h :N h ) 


WREXT, WSEXT( N h , N p ) 

LR, LS( N h , N p ) 

KMATRX( 4*N h *N p , 4*N h *N p ) 

FW( N h ) 

KR, KI( N p , N p ) 

WR, WS( N h , N p ) 

LOAD, WASH( 4*N h *N ) 

LUP, LDN(N h , -N h :N h ) 

LRUP, LRDN, LSUP, LSDN( N h , -N h :N h ) 
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